* Naval Observatory Vector Astronomy Software (NOVAS)
* 
* NOVAS has no licensing requirements.  If you use NOVAS in an
* application, an acknowledgement of the Astronomical Applications
* Department of the U.S. Naval Observatory would be appropriate. Your
* input helps us justify continued development of NOVAS. 
*
* The User's Guide is the official reference for NOVAS C3.1 and may be cited as:
*    Bangert, J., Puatua, W., Kaplan, G., Bartlett, J., Harris, W., Fredericks, A., & Monet, A. 
*    2011, User's Guide to NOVAS Version C3.1 (Washington, DC: USNO).

*   MEMBER 'VASUP1' FOLLOWS
*
      SUBROUTINE APSTAR (TJD,N,RAM,DECM,PMRA,PMDEC,PARLAX,RADVEL,RA,DEC)
      save
C
C     THIS SUBROUTINE COMPUTES THE APPARENT PLACE OF A STAR,
C     GIVEN ITS MEAN PLACE, PROPER MOTION, PARALLAX, AND RADIAL
C     VELOCITY FOR J2000.0.
C
C          TJD    = TDT JULIAN DATE FOR APPARENT PLACE (IN)
C          N      = BODY IDENTIFICATION NUMBER FOR THE EARTH (IN)
C          RAM    = MEAN RIGHT ASCENSION J2000.0 IN HOURS (IN)
C          DECM   = MEAN DECLINATION J2000.0 IN DEGREES (IN)
C          PMRA   = PROPER MOTION IN RA IN SECONDS OF TIME PER JULIAN
C                   CENTURY (IN)
C          PMDEC  = PROPER MOTION IN DEC IN SECONDS OF ARC PER JULIAN
C                   CENTURY (IN)
C          PARLAX = PARALLAX IN SECONDS OF ARC (IN)
C          RADVEL = RADIAL VELOCITY IN KILOMETERS PER SECOND (IN)
C          RA     = APPARENT RIGHT ASCENSION IN HOURS, REFERRED TO
C                   TRUE EQUATOR AND EQUINOX OF DATE (OUT)
C          DEC    = APPARENT DECLINATION IN DEGREES, REFERRED TO
C                   TRUE EQUATOR AND EQUINOX OF DATE (OUT)
C
C
      DOUBLE PRECISION TJD,UJD,RAM,DECM,PMRA,PMDEC,PARLAX,RADVEL,
     /     GLON,GLAT,HT,RA,DEC,T0,T1,TLAST,
     /     X,SECDIF,EQEQ,ST,GAST,RM,DM,PMR,PMD,PI,RV,TLIGHT,R,D,
     /     PEB,VEB,PES,VES,POG,VOG,PB,VB,PS,VS,
     /     POS1,VEL1,POS2,VEL2,POS3,POS4,POS5,POS6,POS7,
     /     DABS,DMOD
      DIMENSION PEB(3), VEB(3), PES(3), VES(3), POG(3), VOG(3),
     /     PB(3), VB(3), PS(3), VS(3),
     /     POS1(3), VEL1(3), POS2(3), VEL2(3),
     /     POS3(3), POS4(3), POS5(3), POS6(3), POS7(3)
C
      DATA T0/2451545.0D0/
C     T0 = TDB JULIAN DATE OF EPOCH J2000.0
      DATA TLAST/0.0D0/
C
C     COMPUTE T1, THE TDB JULIAN DATE CORRESPONDING TO TJD
      CALL TIMES (TJD,X,SECDIF)
      T1 = TJD + SECDIF / 86400.0D0
      IF (DABS(TJD-TLAST).LT.1.0D-6) GO TO 20
C
C     GET POSITION AND VELOCITY OF THE EARTH WRT BARYCENTER OF
C     SOLAR SYSTEM AND WRT CENTER OF SUN
      CALL SOLSYS (T1,N,0,PEB,VEB,IERR)
      IF (IERR.NE.0) GO TO 40
      CALL SOLSYS (T1,N,1,PES,VES,IERR)
      IF (IERR.NE.0) GO TO 40
      TLAST = TJD
C
   20 DO 22 J=1,3
      PB(J) = PEB(J)
      VB(J) = VEB(J)
      PS(J) = PES(J)
      VS(J) = VES(J)
   22 CONTINUE
      RM = RAM
      DM = DECM
      PMR = PMRA
      PMD = PMDEC
      PI = PARLAX
      RV = RADVEL
C
C     COMPUTE APPARENT PLACE
   30 CALL VECTRS (RM,DM,PMR,PMD,PI,RV,POS1,VEL1)
      CALL PROPMO (T0,POS1,VEL1,T1,POS2)
      CALL GEOCEN (POS2,PB,POS3,TLIGHT)
      CALL SUNFLD (POS3,PS,POS4)
      CALL ABERAT (POS4,VB,TLIGHT,POS5)
      CALL PRECES (T0,POS5,T1,POS6)
      CALL NUTATE (T1,POS6,POS7)
      CALL ANGLES (POS7,R,D)
C
      RA = R
      DEC = D
      RETURN
C
   40 RA = 0.0D0
      DEC = 0.0D0
      TLAST = 0.0D0
      RETURN
C
C
      ENTRY TPSTAR (UJD,GLON,GLAT,HT,RA,DEC)
C
C     THIS ENTRY COMPUTES THE TOPOCENTRIC PLACE OF A STAR,
C     GIVEN THE LOCATION OF THE OBSERVER.  THIS ENTRY ASSUMES APSTAR
C     WAS PREVIOUSLY CALLED, AND USES DATA COMPUTED BY APSTAR.
C
C          UJD    = UT1 JULIAN DATE, OR EQUIVALENT GREENWICH APPARENT
C                   SIDEREAL TIME IN HOURS, FOR TOPOCENTRIC PLACE (IN)
C          GLON   = GEODETIC LONGITUDE (EAST +) OF OBSERVER
C                   IN DEGREES (IN)
C          GLAT   = GEODETIC LATITUDE (NORTH +) OF OBSERVER
C                   IN DEGREES (IN)
C          HT     = HEIGHT OF OBSERVER IN METERS (IN)
C          RA     = TOPOCENTRIC RIGHT ASCENSION IN HOURS, REFERRED TO
C                   TRUE EQUATOR AND EQUINOX OF DATE (OUT)
C          DEC    = TOPOCENTRIC DECLINATION IN DEGREES, REFERRED TO
C                   TRUE EQUATOR AND EQUINOX OF DATE (OUT)
C
C
C     GET POSITION AND VELOCITY OF OBSERVER WRT CENTER OF EARTH
   50 IF (UJD.GT.100.0D0) GO TO 52
      GAST = DMOD(UJD,24.0D0)
      GO TO 55
   52 CALL SIDTIM (UJD,0.0D0,0,ST)
      CALL ETILT (T1,X,X,EQEQ,X,X)
      GAST = ST + EQEQ/3600.0D0
   55 CALL TERRA (GLON,GLAT,HT,GAST,POS1,VEL1)
      CALL NUTATE (-T1,POS1,POS2)
      CALL PRECES (T1,POS2,T0,POG)
      CALL NUTATE (-T1,VEL1,VEL2)
      CALL PRECES (T1,VEL2,T0,VOG)
C
C     COMPUTE POSITION AND VELOCITY OF OBSERVER WRT BARYCENTER OF
C     SOLAR SYSTEM AND WRT CENTER OF SUN
   60 DO 62 J=1,3
      PB(J) = PEB(J) + POG(J)
      VB(J) = VEB(J) + VOG(J)
      PS(J) = PES(J) + POG(J)
      VS(J) = VES(J) + VOG(J)
   62 CONTINUE
C
C     RECOMPUTE APPARENT PLACE USING POSITION AND VELOCITY OF OBSERVER
      GO TO 30
C
      END
      SUBROUTINE APPLAN (TJD,L,N,RA,DEC,DIS)
      save
C
C     THIS SUBROUTINE COMPUTES THE APPARENT PLACE OF A PLANET OR OTHER
C     SOLAR SYSTEM BODY.  RECTANGULAR COORDINATES OF SOLAR SYSTEM BODIES
C     ARE OBTAINED FROM SUBROUTINE SOLSYS.
C
C          TJD    = TDT JULIAN DATE FOR APPARENT PLACE (IN)
C          L      = BODY IDENTIFICATION NUMBER FOR DESIRED PLANET (IN)
C          N      = BODY IDENTIFICATION NUMBER FOR THE EARTH (IN)
C          RA     = APPARENT RIGHT ASCENSION IN HOURS, REFERRED TO
C                   TRUE EQUATOR AND EQUINOX OF DATE (OUT)
C          DEC    = APPARENT DECLINATION IN DEGREES, REFERRED TO
C                   TRUE EQUATOR AND EQUINOX OF DATE (OUT)
C          DIS    = TRUE DISTANCE FROM EARTH TO PLANET IN AU (OUT)
C
C
      DOUBLE PRECISION TJD,UJD,GLON,GLAT,HT,RA,DEC,DIS,
     /     T0,T1,T2,T3,TLAST,C,
     /     X,SECDIF,EQEQ,ST,GAST,TLIGHT,R,D,S,
     /     PEB,VEB,PES,VES,POG,VOG,PB,VB,PS,VS,
     /     POS1,VEL1,POS2,VEL2,POS3,POS4,POS5,POS6,
     /     DABS,DMOD
      DIMENSION PEB(3), VEB(3), PES(3), VES(3), POG(3), VOG(3),
     /     PB(3), VB(3), PS(3), VS(3),
     /     POS1(3), VEL1(3), POS2(3), VEL2(3),
     /     POS3(3), POS4(3), POS5(3), POS6(3)
C
      DATA C/173.1446333D0/
C     C = SPEED OF LIGHT IN AU/DAY = 86400 / 499.004782
      DATA T0/2451545.0D0/
C     T0 = TDB JULIAN DATE OF EPOCH J2000.0
      DATA TLAST/0.0D0/
C
C     COMPUTE T1, THE TDB JULIAN DATE CORRESPONDING TO TJD
      CALL TIMES (TJD,X,SECDIF)
      T1 = TJD + SECDIF / 86400.0D0
      IF (DABS(TJD-TLAST).LT.1.0D-6) GO TO 20
C
C     GET POSITION AND VELOCITY OF THE EARTH WRT BARYCENTER OF
C     SOLAR SYSTEM AND WRT CENTER OF SUN
      CALL SOLSYS (T1,N,0,PEB,VEB,IERR)
      IF (IERR.NE.0) GO TO 40
      CALL SOLSYS (T1,N,1,PES,VES,IERR)
      IF (IERR.NE.0) GO TO 40
      TLAST = TJD
C
   20 DO 22 J=1,3
      PB(J) = PEB(J)
      VB(J) = VEB(J)
      PS(J) = PES(J)
      VS(J) = VES(J)
   22 CONTINUE
      LPLAN = L
C
C     GET POSITION OF PLANET WRT BARYCENTER OF SOLAR SYSTEM
   30 CALL SOLSYS (T1,LPLAN,0,POS1,VEL1,IERR)
      IF (IERR.NE.0) GO TO 40
      CALL GEOCEN (POS1,PB,POS2,TLIGHT)
      S = TLIGHT * C
      T2 = T1 - TLIGHT
   33 CALL SOLSYS (T2,LPLAN,0,POS1,VEL1,IERR)
      IF (IERR.NE.0) GO TO 40
      CALL GEOCEN (POS1,PB,POS2,TLIGHT)
      T3 = T1 - TLIGHT
      IF (DABS(T3-T2).LT.1.0D-8) GO TO 35
      T2 = T3
      GO TO 33
C
C     FINISH APPARENT PLACE COMPUTATION
   35 CONTINUE
      CALL SUNFLD (POS2,PS,POS3)
      CALL ABERAT (POS3,VB,TLIGHT,POS4)
      CALL PRECES (T0,POS4,T1,POS5)
      CALL NUTATE (T1,POS5,POS6)
      CALL ANGLES (POS6,R,D)
      RA = R
      DEC = D
      DIS = S
      RETURN
C
   40 RA = 0.0D0
      DEC = 0.0D0
      DIS = 0.0D0
      TLAST = 0.0D0
      RETURN
C
C
      ENTRY TPPLAN (UJD,GLON,GLAT,HT,RA,DEC,DIS)
C
C     THIS ENTRY COMPUTES THE TOPOCENTRIC PLACE OF A PLANET,
C     GIVEN THE LOCATION OF THE OBSERVER.  THIS ENTRY ASSUMES APPLAN
C     WAS PREVIOUSLY CALLED, AND USES DATA COMPUTED BY APPLAN.
C
C          UJD    = UT1 JULIAN DATE, OR EQUIVALENT GREENWICH APPARENT
C                   SIDEREAL TIME IN HOURS, FOR TOPOCENTRIC PLACE (IN)
C          GLON   = GEODETIC LONGITUDE (EAST +) OF OBSERVER
C                   IN DEGREES (IN)
C          GLAT   = GEODETIC LATITUDE (NORTH +) OF OBSERVER
C                   IN DEGREES (IN)
C          HT     = HEIGHT OF OBSERVER IN METERS (IN)
C          RA     = TOPOCENTRIC RIGHT ASCENSION IN HOURS, REFERRED TO
C                   TRUE EQUATOR AND EQUINOX OF DATE (OUT)
C          DEC    = TOPOCENTRIC DECLINATION IN DEGREES, REFERRED TO
C                   TRUE EQUATOR AND EQUINOX OF DATE (OUT)
C          DIS    = TRUE DISTANCE FROM OBSERVER TO PLANET IN AU (OUT)
C
C
C     GET POSITION AND VELOCITY OF OBSERVER WRT CENTER OF EARTH
   50 IF (UJD.GT.100.0D0) GO TO 52
      GAST = DMOD(UJD,24.0D0)
      GO TO 55
   52 CALL SIDTIM (UJD,0.0D0,0,ST)
      CALL ETILT (T1,X,X,EQEQ,X,X)
      GAST = ST + EQEQ/3600.0D0
   55 CALL TERRA (GLON,GLAT,HT,GAST,POS1,VEL1)
      CALL NUTATE (-T1,POS1,POS2)
      CALL PRECES (T1,POS2,T0,POG)
      CALL NUTATE (-T1,VEL1,VEL2)
      CALL PRECES (T1,VEL2,T0,VOG)
C
C     COMPUTE POSITION AND VELOCITY OF OBSERVER WRT BARYCENTER OF
C     SOLAR SYSTEM AND WRT CENTER OF SUN
   60 DO 62 J=1,3
      PB(J) = PEB(J) + POG(J)
      VB(J) = VEB(J) + VOG(J)
      PS(J) = PES(J) + POG(J)
      VS(J) = VES(J) + VOG(J)
   62 CONTINUE
C
C     RECOMPUTE APPARENT PLACE USING POSITION AND VELOCITY OF OBSERVER
      GO TO 30
C
      END
      SUBROUTINE VPSTAR (TJD,N,RAM,DECM,PMRA,PMDEC,PARLAX,RADVEL,RA,DEC)
      save
C
C     THIS SUBROUTINE COMPUTES THE VIRTUAL PLACE OF A STAR,
C     GIVEN ITS MEAN PLACE, PROPER MOTION, PARALLAX, AND RADIAL
C     VELOCITY FOR J2000.0.
C
C          TJD    = TDT JULIAN DATE FOR VIRTUAL PLACE (IN)
C          N      = BODY IDENTIFICATION NUMBER FOR THE EARTH (IN)
C          RAM    = MEAN RIGHT ASCENSION J2000.0 IN HOURS (IN)
C          DECM   = MEAN DECLINATION J2000.0 IN DEGREES (IN)
C          PMRA   = PROPER MOTION IN RA IN SECONDS OF TIME PER JULIAN
C                   CENTURY (IN)
C          PMDEC  = PROPER MOTION IN DEC IN SECONDS OF ARC PER JULIAN
C                   CENTURY (IN)
C          PARLAX = PARALLAX IN SECONDS OF ARC (IN)
C          RADVEL = RADIAL VELOCITY IN KILOMETERS PER SECOND (IN)
C          RA     = VIRTUAL RIGHT ASCENSION IN HOURS, REFERRED TO
C                   MEAN EQUATOR AND EQUINOX OF J2000.0 (OUT)
C          DEC    = VIRTUAL DECLINATION IN DEGREES, REFERRED TO
C                   MEAN EQUATOR AND EQUINOX OF J2000.0 (OUT)
C
C
      DOUBLE PRECISION TJD,UJD,RAM,DECM,PMRA,PMDEC,PARLAX,RADVEL,
     /     GLON,GLAT,HT,RA,DEC,T0,T1,TLAST,
     /     X,SECDIF,EQEQ,ST,GAST,RM,DM,PMR,PMD,PI,RV,TLIGHT,R,D,
     /     PEB,VEB,PES,VES,POG,VOG,PB,VB,PS,VS,
     /     POS1,VEL1,POS2,VEL2,POS3,POS4,POS5,
     /     DABS,DMOD
      DIMENSION PEB(3), VEB(3), PES(3), VES(3), POG(3), VOG(3),
     /     PB(3), VB(3), PS(3), VS(3),
     /     POS1(3), VEL1(3), POS2(3), VEL2(3),
     /     POS3(3), POS4(3), POS5(3)
C
      DATA T0/2451545.0D0/
C     T0 = TDB JULIAN DATE OF EPOCH J2000.0
      DATA TLAST/0.0D0/
C
C     COMPUTE T1, THE TDB JULIAN DATE CORRESPONDING TO TJD
      CALL TIMES (TJD,X,SECDIF)
      T1 = TJD + SECDIF / 86400.0D0
      IF (DABS(TJD-TLAST).LT.1.0D-6) GO TO 20
C
C     GET POSITION AND VELOCITY OF THE EARTH WRT BARYCENTER OF
C     SOLAR SYSTEM AND WRT CENTER OF SUN
      CALL SOLSYS (T1,N,0,PEB,VEB,IERR)
      IF (IERR.NE.0) GO TO 40
      CALL SOLSYS (T1,N,1,PES,VES,IERR)
      IF (IERR.NE.0) GO TO 40
      TLAST = TJD
C
   20 DO 22 J=1,3
      PB(J) = PEB(J)
      VB(J) = VEB(J)
      PS(J) = PES(J)
      VS(J) = VES(J)
   22 CONTINUE
      RM = RAM
      DM = DECM
      PMR = PMRA
      PMD = PMDEC
      PI = PARLAX
      RV = RADVEL
C
C     COMPUTE VIRTUAL PLACE
   30 CALL VECTRS (RM,DM,PMR,PMD,PI,RV,POS1,VEL1)
      CALL PROPMO (T0,POS1,VEL1,T1,POS2)
      CALL GEOCEN (POS2,PB,POS3,TLIGHT)
      CALL SUNFLD (POS3,PS,POS4)
      CALL ABERAT (POS4,VB,TLIGHT,POS5)
      CALL ANGLES (POS5,R,D)
C
      RA = R
      DEC = D
      RETURN
C
   40 RA = 0.0D0
      DEC = 0.0D0
      TLAST = 0.0D0
      RETURN
C
C
      ENTRY LPSTAR (UJD,GLON,GLAT,HT,RA,DEC)
C
C     THIS ENTRY COMPUTES THE LOCAL PLACE OF A STAR,
C     GIVEN THE LOCATION OF THE OBSERVER.  THIS ENTRY ASSUMES VPSTAR
C     WAS PREVIOUSLY CALLED, AND USES DATA COMPUTED BY VPSTAR.
C
C          UJD    = UT1 JULIAN DATE, OR EQUIVALENT GREENWICH APPARENT
C                   SIDEREAL TIME IN HOURS, FOR LOCAL PLACE (IN)
C          GLON   = GEODETIC LONGITUDE (EAST +) OF OBSERVER
C                   IN DEGREES (IN)
C          GLAT   = GEODETIC LATITUDE (NORTH +) OF OBSERVER
C                   IN DEGREES (IN)
C          HT     = HEIGHT OF OBSERVER IN METERS (IN)
C          RA     = LOCAL RIGHT ASCENSION IN HOURS, REFERRED TO
C                   MEAN EQUATOR AND EQUINOX OF J2000.0 (OUT)
C          DEC    = LOCAL DECLINATION IN DEGREES, REFERRED TO
C                   MEAN EQUATOR AND EQUINOX OF J2000.0 (OUT)
C
C
C     GET POSITION AND VELOCITY OF OBSERVER WRT CENTER OF EARTH
   50 IF (UJD.GT.100.0D0) GO TO 52
      GAST = DMOD(UJD,24.0D0)
      GO TO 55
   52 CALL SIDTIM (UJD,0.0D0,0,ST)
      CALL ETILT (T1,X,X,EQEQ,X,X)
      GAST = ST + EQEQ/3600.0D0
   55 CALL TERRA (GLON,GLAT,HT,GAST,POS1,VEL1)
      CALL NUTATE (-T1,POS1,POS2)
      CALL PRECES (T1,POS2,T0,POG)
      CALL NUTATE (-T1,VEL1,VEL2)
      CALL PRECES (T1,VEL2,T0,VOG)
C
C     COMPUTE POSITION AND VELOCITY OF OBSERVER WRT BARYCENTER OF
C     SOLAR SYSTEM AND WRT CENTER OF SUN
   60 DO 62 J=1,3
      PB(J) = PEB(J) + POG(J)
      VB(J) = VEB(J) + VOG(J)
      PS(J) = PES(J) + POG(J)
      VS(J) = VES(J) + VOG(J)
   62 CONTINUE
C
C     RECOMPUTE VIRTUAL PLACE USING POSITION AND VELOCITY OF OBSERVER
      GO TO 30
C
      END
      SUBROUTINE VPPLAN (TJD,L,N,RA,DEC,DIS)
      save
C
C     THIS SUBROUTINE COMPUTES THE VIRTUAL PLACE OF A PLANET OR OTHER
C     SOLAR SYSTEM BODY.  RECTANGULAR COORDINATES OF SOLAR SYSTEM BODIES
C     ARE OBTAINED FROM SUBROUTINE SOLSYS.
C
C          TJD    = TDT JULIAN DATE FOR VIRTUAL PLACE (IN)
C          L      = BODY IDENTIFICATION NUMBER FOR DESIRED PLANET (IN)
C          N      = BODY IDENTIFICATION NUMBER FOR THE EARTH (IN)
C          RA     = VIRTUAL RIGHT ASCENSION IN HOURS, REFERRED TO
C                   MEAN EQUATOR AND EQUINOX OF J2000.0 (OUT)
C          DEC    = VIRTUAL DECLINATION IN DEGREES, REFERRED TO
C                   MEAN EQUATOR AND EQUINOX OF J2000.0 (OUT)
C          DIS    = TRUE DISTANCE FROM EARTH TO PLANET IN AU (OUT)
C
C
      DOUBLE PRECISION TJD,UJD,GLON,GLAT,HT,RA,DEC,DIS,
     /     T0,T1,T2,T3,TLAST,C,
     /     X,SECDIF,EQEQ,ST,GAST,TLIGHT,R,D,S,
     /     PEB,VEB,PES,VES,POG,VOG,PB,VB,PS,VS,
     /     POS1,VEL1,POS2,VEL2,POS3,POS4,
     /     DABS,DMOD
      DIMENSION PEB(3), VEB(3), PES(3), VES(3), POG(3), VOG(3),
     /     PB(3), VB(3), PS(3), VS(3),
     /     POS1(3), VEL1(3), POS2(3), VEL2(3),
     /     POS3(3), POS4(3)
C
      DATA C/173.1446333D0/
C     C = SPEED OF LIGHT IN AU/DAY = 86400 / 499.004782
      DATA T0/2451545.0D0/
C     T0 = TDB JULIAN DATE OF EPOCH J2000.0
      DATA TLAST/0.0D0/
C
C     COMPUTE T1, THE TDB JULIAN DATE CORRESPONDING TO TJD
      CALL TIMES (TJD,X,SECDIF)
      T1 = TJD + SECDIF / 86400.0D0
      IF (DABS(TJD-TLAST).LT.1.0D-6) GO TO 20
C
C     GET POSITION AND VELOCITY OF THE EARTH WRT BARYCENTER OF
C     SOLAR SYSTEM AND WRT CENTER OF SUN
      CALL SOLSYS (T1,N,0,PEB,VEB,IERR)
      IF (IERR.NE.0) GO TO 40
      CALL SOLSYS (T1,N,1,PES,VES,IERR)
      IF (IERR.NE.0) GO TO 40
      TLAST = TJD
C
   20 DO 22 J=1,3
      PB(J) = PEB(J)
      VB(J) = VEB(J)
      PS(J) = PES(J)
      VS(J) = VES(J)
   22 CONTINUE
      LPLAN = L
C
C     GET POSITION OF PLANET WRT BARYCENTER OF SOLAR SYSTEM
   30 CALL SOLSYS (T1,LPLAN,0,POS1,VEL1,IERR)
      IF (IERR.NE.0) GO TO 40
      CALL GEOCEN (POS1,PB,POS2,TLIGHT)
      S = TLIGHT * C
      T2 = T1 - TLIGHT
   33 CALL SOLSYS (T2,LPLAN,0,POS1,VEL1,IERR)
      IF (IERR.NE.0) GO TO 40
      CALL GEOCEN (POS1,PB,POS2,TLIGHT)
      T3 = T1 - TLIGHT
      IF (DABS(T3-T2).LT.1.0D-8) GO TO 35
      T2 = T3
      GO TO 33
C
C     FINISH VIRTUAL PLACE COMPUTATION
   35 CONTINUE
      CALL SUNFLD (POS2,PS,POS3)
      CALL ABERAT (POS3,VB,TLIGHT,POS4)
      CALL ANGLES (POS4,R,D)
      RA = R
      DEC = D
      DIS = S
      RETURN
C
   40 RA = 0.0D0
      DEC = 0.0D0
      DIS = 0.0D0
      TLAST = 0.0D0
      RETURN
C
C
      ENTRY LPPLAN (UJD,GLON,GLAT,HT,RA,DEC,DIS)
C
C     THIS ENTRY COMPUTES THE LOCAL PLACE OF A PLANET, GIVEN
C     THE LOCATION OF THE OBSERVER.  THIS ENTRY ASSUMES VPPLAN WAS
C     PREVIOUSLY CALLED, AND USES DATA COMPUTED BY VPPLAN.
C
C          UJD    = UT1 JULIAN DATE, OR EQUIVALENT GREENWICH APPARENT
C                   SIDEREAL TIME IN HOURS, FOR LOCAL PLACE (IN)
C          GLON   = GEODETIC LONGITUDE (EAST +) OF OBSERVER
C                   IN DEGREES (IN)
C          GLAT   = GEODETIC LATITUDE (NORTH +) OF OBSERVER
C                   IN DEGREES (IN)
C          HT     = HEIGHT OF OBSERVER IN METERS (IN)
C          RA     = LOCAL RIGHT ASCENSION IN HOURS, REFERRED TO
C                   MEAN EQUATOR AND EQUINOX OF J2000.0 (OUT)
C          DEC    = LOCAL DECLINATION IN DEGREES, REFERRED TO
C                   MEAN EQUATOR AND EQUINOX OF J2000.0 (OUT)
C          DIS    = TRUE DISTANCE FROM OBSERVER TO PLANET IN AU (OUT)
C
C
C     GET POSITION AND VELOCITY OF OBSERVER WRT CENTER OF EARTH
   50 IF (UJD.GT.100.0D0) GO TO 52
      GAST = DMOD(UJD,24.0D0)
      GO TO 55
   52 CALL SIDTIM (UJD,0.0D0,0,ST)
      CALL ETILT (T1,X,X,EQEQ,X,X)
      GAST = ST + EQEQ/3600.0D0
   55 CALL TERRA (GLON,GLAT,HT,GAST,POS1,VEL1)
      CALL NUTATE (-T1,POS1,POS2)
      CALL PRECES (T1,POS2,T0,POG)
      CALL NUTATE (-T1,VEL1,VEL2)
      CALL PRECES (T1,VEL2,T0,VOG)
C
C     COMPUTE POSITION AND VELOCITY OF OBSERVER WRT BARYCENTER OF
C     SOLAR SYSTEM AND WRT CENTER OF SUN
   60 DO 62 J=1,3
      PB(J) = PEB(J) + POG(J)
      VB(J) = VEB(J) + VOG(J)
      PS(J) = PES(J) + POG(J)
      VS(J) = VES(J) + VOG(J)
   62 CONTINUE
C
C     RECOMPUTE VIRTUAL PLACE USING POSITION AND VELOCITY OF OBSERVER
      GO TO 30
C
      END
      SUBROUTINE ASSTAR (TJD,N,RAM,DECM,PMRA,PMDEC,PARLAX,RADVEL,RA,DEC)
      save
C
C     THIS SUBROUTINE COMPUTES THE ASTROMETRIC PLACE OF A STAR,
C     GIVEN ITS MEAN PLACE, PROPER MOTION, PARALLAX, AND RADIAL
C     VELOCITY FOR J2000.0.
C
C          TJD    = TDT JULIAN DATE FOR ASTROMETRIC PLACE (IN)
C          N      = BODY IDENTIFICATION NUMBER FOR THE EARTH (IN)
C          RAM    = MEAN RIGHT ASCENSION J2000.0 IN HOURS (IN)
C          DECM   = MEAN DECLINATION J2000.0 IN DEGREES (IN)
C          PMRA   = PROPER MOTION IN RA IN SECONDS OF TIME PER JULIAN
C                   CENTURY (IN)
C          PMDEC  = PROPER MOTION IN DEC IN SECONDS OF ARC PER JULIAN
C                   CENTURY (IN)
C          PARLAX = PARALLAX IN SECONDS OF ARC (IN)
C          RADVEL = RADIAL VELOCITY IN KILOMETERS PER SECOND (IN)
C          RA     = ASTROMETRIC RIGHT ASCENSION IN HOURS, REFERRED TO
C                   MEAN EQUATOR AND EQUINOX OF J2000.0 (OUT)
C          DEC    = ASTROMETRIC DECLINATION IN DEGREES, REFERRED TO
C                   MEAN EQUATOR AND EQUINOX OF J2000.0 (OUT)
C
C
      DOUBLE PRECISION TJD,RAM,DECM,PMRA,PMDEC,PARLAX,RADVEL,RA,DEC,
     /     T0,T1,TLAST,X,SECDIF,RM,DM,PMR,PMD,PI,RV,TLIGHT,R,D,
     /     PEB,VEB,PB,VB,POS1,VEL1,POS2,POS3,DABS
      DIMENSION PEB(3), VEB(3), PB(3), VB(3),
     /     POS1(3), VEL1(3), POS2(3), POS3(3)
C
      DATA T0/2451545.0D0/
C     T0 = TDB JULIAN DATE OF EPOCH J2000.0
      DATA TLAST/0.0D0/
C
C     COMPUTE T1, THE TDB JULIAN DATE CORRESPONDING TO TJD
      CALL TIMES (TJD,X,SECDIF)
      T1 = TJD + SECDIF / 86400.0D0
      IF (DABS(TJD-TLAST).LT.1.0D-6) GO TO 20
C
C     GET POSITION AND VELOCITY OF THE EARTH WRT BARYCENTER OF
C     SOLAR SYSTEM
      CALL SOLSYS (T1,N,0,PEB,VEB,IERR)
      IF (IERR.NE.0) GO TO 40
      TLAST = TJD
C
   20 DO 22 J=1,3
      PB(J) = PEB(J)
      VB(J) = VEB(J)
   22 CONTINUE
      RM = RAM
      DM = DECM
      PMR = PMRA
      PMD = PMDEC
      PI = PARLAX
      RV = RADVEL
C
C     COMPUTE ASTROMETRIC PLACE
   30 CALL VECTRS (RM,DM,PMR,PMD,PI,RV,POS1,VEL1)
      CALL PROPMO (T0,POS1,VEL1,T1,POS2)
      CALL GEOCEN (POS2,PB,POS3,TLIGHT)
      CALL ANGLES (POS3,R,D)
      RA = R
      DEC = D
      RETURN
C
   40 RA = 0.0D0
      DEC = 0.0D0
      TLAST = 0.0D0
      RETURN
      END
      SUBROUTINE ASPLAN (TJD,L,N,RA,DEC,DIS)
      save
C
C     THIS SUBROUTINE COMPUTES THE ASTROMETRIC PLACE OF A PLANET OR
C     OTHER SOLAR SYSTEM BODY.   RECTANGULAR COORDINATES OF SOLAR SYSTEM
C     BODIES ARE OBTAINED FROM SUBROUTINE SOLSYS.
C
C          TJD    = TDT JULIAN DATE FOR ASTROMETRIC PLACE (IN)
C          L      = BODY IDENTIFICATION NUMBER FOR DESIRED PLANET (IN)
C          N      = BODY IDENTIFICATION NUMBER FOR THE EARTH (IN)
C          RA     = ASTROMETRIC RIGHT ASCENSION IN HOURS, REFERRED TO
C                   MEAN EQUATOR AND EQUINOX OF J2000.0 (OUT)
C          DEC    = ASTROMETRIC DECLINATION IN DEGREES, REFERRED TO
C                   MEAN EQUATOR AND EQUINOX OF J2000.0 (OUT)
C          DIS    = TRUE DISTANCE FROM EARTH TO PLANET IN AU (OUT)
C
C
      DOUBLE PRECISION TJD,RA,DEC,DIS,T0,T1,T2,T3,TLAST,X,SECDIF,
     /     C,TLIGHT,R,D,S,
     /     PEB,VEB,PB,VB,POS1,VEL1,POS2,DABS
      DIMENSION PEB(3), VEB(3), PB(3), VB(3),
     /     POS1(3), VEL1(3), POS2(3)
C
      DATA C/173.1446333D0/
C     C = SPEED OF LIGHT IN AU/DAY = 86400 / 499.004782
      DATA T0/2451545.0D0/
C     T0 = TDB JULIAN DATE OF EPOCH J2000.0
      DATA TLAST/0.0D0/
C
C     COMPUTE T1, THE TDB JULIAN DATE CORRESPONDING TO TJD
      CALL TIMES (TJD,X,SECDIF)
      T1 = TJD + SECDIF / 86400.0D0
      IF (DABS(TJD-TLAST).LT.1.0D-6) GO TO 20
C
C     GET POSITION AND VELOCITY OF THE EARTH WRT BARYCENTER OF
C     SOLAR SYSTEM
      CALL SOLSYS (T1,N,0,PEB,VEB,IERR)
      IF (IERR.NE.0) GO TO 40
      TLAST = TJD
C
   20 DO 22 J=1,3
      PB(J) = PEB(J)
      VB(J) = VEB(J)
   22 CONTINUE
      LPLAN = L
C
C     GET POSITION OF PLANET WRT BARYCENTER OF SOLAR SYSTEM
   30 CALL SOLSYS (T1,LPLAN,0,POS1,VEL1,IERR)
      IF (IERR.NE.0) GO TO 40
      CALL GEOCEN (POS1,PB,POS2,TLIGHT)
      S = TLIGHT * C
      T2 = T1 - TLIGHT
   33 CALL SOLSYS (T2,LPLAN,0,POS1,VEL1,IERR)
      IF (IERR.NE.0) GO TO 40
      CALL GEOCEN (POS1,PB,POS2,TLIGHT)
      T3 = T1 - TLIGHT
      IF (DABS(T3-T2).LT.1.0D-8) GO TO 35
      T2 = T3
      GO TO 33
C
C     FINISH ASTROMETRIC PLACE COMPUTATION
   35 CONTINUE
      CALL ANGLES (POS2,R,D)
      RA = R
      DEC = D
      DIS = S
      RETURN
C
   40 RA = 0.0D0
      DEC = 0.0D0
      DIS = 0.0D0
      TLAST = 0.0D0
      RETURN
      END
      SUBROUTINE MPSTAR (TJD,N,RA,DEC,RAM,DECM)
      save
C
C     THIS SUBROUTINE COMPUTES THE MEAN PLACE OF A STAR FOR J2000.0,
C     GIVEN ITS APPARENT PLACE AT DATE TJD.  PROPER MOTION, PARALLAX,
C     AND RADIAL VELOCITY ARE ASSUMED TO BE ZERO.
C
C          TJD    = TDT JULIAN DATE OF APPARENT PLACE (IN)
C          N      = BODY IDENTIFICATION NUMBER FOR THE EARTH (IN)
C          RA     = APPARENT RIGHT ASCENSION IN HOURS, REFERRED TO
C                   TRUE EQUATOR AND EQUINOX OF DATE (IN)
C          DEC    = APPARENT DECLINATION IN DEGREES, REFERRED TO
C                   TRUE EQUATOR AND EQUINOX OF DATE (IN)
C          RAM    = MEAN RIGHT ASCENSION J2000.0 IN HOURS (OUT)
C          DECM   = MEAN DECLINATION J2000.0 IN DEGREES (OUT)
C
C
      DOUBLE PRECISION TJD,RA,DEC,RAM,DECM,T1,RAMNEW,DCMNEW,
     /     RAMOLD,DCMOLD,R,D,DELRA,DELDEC,DABS,DMOD
C
      T1 = TJD
      RAMNEW = DMOD(RA,24.0D0)
      IF (RAMNEW.LT.0.0D0) RAMNEW = RAMNEW + 24.0D0
      DCMNEW = DEC
      ITER = 0
C
   20 ITER = ITER + 1
      RAMOLD = RAMNEW
      DCMOLD = DCMNEW
      R = RAMOLD
      D = DCMOLD
      CALL APSTAR (T1,N,R,D,0.0D0,0.0D0,0.0D0,0.0D0,R,D)
      DELRA = R - RAMOLD
      DELDEC = D - DCMOLD
      IF (DELRA.LT.-12.0D0) DELRA = DELRA + 24.0D0
      IF (DELRA.GT.+12.0D0) DELRA = DELRA - 24.0D0
      RAMNEW = RA - DELRA
      DCMNEW = DEC - DELDEC
      IF (ITER.GT.20) GO TO 40
      IF (DABS(RAMNEW-RAMOLD).GT.1.0D-10) GO TO 20
      IF (DABS(DCMNEW-DCMOLD).GT.1.0D-09) GO TO 20
C
      RAM = RAMNEW
      DECM = DCMNEW
      IF (RAM.LT. 0.0D0) RAM = RAM + 24.0D0
      IF (RAM.GE.24.0D0) RAM = RAM - 24.0D0
      RETURN
C
   40 RAM = 0.0D0
      DECM = 0.0D0
      RETURN
      END
      SUBROUTINE SIDTIM (TJDH,TJDL,K,GST)
      save
C
C     THIS SUBROUTINE COMPUTES THE GREENWICH SIDEREAL TIME
C     (EITHER MEAN OR APPARENT) AT JULIAN DATE TJDH + TJDL.
C     SEE AOKI, ET AL. (1982) ASTRONOMY AND ASTROPYSICS 105, 359-361.
C
C          TJDH   = JULIAN DATE, HIGH-ORDER PART (IN)
C          TJDL   = JULIAN DATE, LOW-ORDER PART (IN)
C                   JULIAN DATE MAY BE SPLIT AT ANY POINT, BUT
C                   FOR HIGHEST PRECISION, SET TJDH TO BE THE INTEGRAL
C                   PART OF THE JULIAN DATE, AND SET TJDL TO BE THE
C                   FRACTIONAL PART
C          K      = TIME SELECTION CODE (IN)
C                   SET K=0 FOR GREENWICH MEAN SIDEREAL TIME
C                   SET K=1 FOR GREENWICH APPARENT SIDEREAL TIME
C          GST    = GREENWICH (MEAN OR APPARENT) SIDEREAL TIME
C                   IN HOURS (OUT)
C
C
      DOUBLE PRECISION TJDH,TJDL,TJD,TH,TL,T0,T,T2,T3,GST,
     /     X,EQEQ,ST,DMOD
C
      DATA T0/2451545.0D0/
C     T0 = TDB JULIAN DATE OF EPOCH J2000.0
C
      TJD = TJDH + TJDL
      TH = (TJDH - T0) / 36525.0D0
      TL =  TJDL       / 36525.0D0
      T = TH + TL
      T2 = T * T
      T3 = T2 * T
C
C     FOR APPARENT SIDEREAL TIME, OBTAIN EQUATION OF THE EQUINOXES
      EQEQ = 0.0D0
      IF (K.EQ.1) CALL ETILT (TJD,X,X,EQEQ,X,X)
C
      ST = EQEQ - 6.2D-6*T3 + 0.093104D0*T2 + 67310.54841D0
     /     + 8640184.812866D0 *TL
     /     + 3155760000.0D0   *TL
     /     + 8640184.812866D0 *TH
     /     + 3155760000.0D0   *TH
C
      GST = DMOD (ST / 3600.0D0, 24.0D0)
      IF (GST.LT.0.0D0) GST = GST + 24.0D0
      RETURN
      END
      SUBROUTINE PNSW (TJD,GAST,X,Y,VECE,VECS)
      save
C
C     TRANSFORMS A VECTOR FROM EARTH-FIXED SYSTEM TO SPACE-FIXED SYSTEM
C     BY APPLYING ROTATIONS FOR WOBBLE, SPIN, NUTATION, AND PRECESSION.
C     (COMBINED ROTATION IS SYMBOLIZED  P N S W .)   SPECIFICALLY,
C     IT TRANSFORMS A VECTOR FROM EARTH-FIXED GEOGRAPHIC SYSTEM TO
C     SPACE-FIXED SYSTEM BASED ON MEAN EQUATOR AND EQUINOX OF J2000.0.
C
C          TJD    = TDT JULIAN DATE (IN)
C          GAST   = GREENWICH APPARENT SIDEREAL TIME, IN HOURS (IN)
C          X      = CONVENTIONALLY-DEFINED X COORDINATE OF ROTATIONAL
C                   POLE WITH RESPECT TO CIO, IN ARCSECONDS (IN)
C          Y      = CONVENTIONALLY-DEFINED Y COORDINATE OF ROTATIONAL
C                   POLE WITH RESPECT TO CIO, IN ARCSECONDS (IN)
C          VECE   = VECTOR IN GEOCENTRIC RECTANGULAR
C                   EARTH-FIXED SYSTEM, REFERRED TO GEOGRAPHIC
C                   EQUATOR AND GREENWICH MERIDIAN (IN)
C          VECS   = VECTOR IN GEOCENTRIC RECTANGULAR
C                   SPACE-FIXED SYSTEM, REFERRED TO MEAN EQUATOR
C                   AND EQUINOX OF J2000.0 (OUT)
C
C          NOTE --  TJD=0 MEANS NO PRECESSION/NUTATION TRANSFORMATION
C                   GAST=0 MEANS NO SPIN TRANSFORMATION
C                   X=Y=0 MEANS NO WOBBLE TRANSFORMATION
C
C
      DOUBLE PRECISION TJD,GAST,X,Y,VECE,VECS,T0,T1,Z,SECDIF,
     /     V1,V2,V3
      DIMENSION VECE(3), VECS(3), V1(3), V2(3), V3(3)
C
      DATA T0/2451545.0D0/
C     T0 = TDB JULIAN DATE OF EPOCH J2000.0
C
C     COMPUTE T1, THE TDB JULIAN DATE CORRESPONDING TO TJD
      IF (TJD.EQ.0.0D0) GO TO 20
      CALL TIMES (TJD,Z,SECDIF)
      T1 = TJD + SECDIF / 86400.0D0
C
   20 IF (X.EQ.0.0D0 .AND. Y.EQ.0.0D0) GO TO 25
      CALL WOBBLE (X,Y,VECE,V1)
      GO TO 30
   25 DO 28 J=1,3
   28 V1(J) = VECE(J)
C
   30 IF (GAST.EQ.0.0D0) GO TO 35
      CALL SPIN (GAST,V1,V2)
      GO TO 40
   35 DO 38 J=1,3
   38 V2(J) = V1(J)
C
   40 IF (TJD.EQ.0.0D0) GO TO 45
      CALL NUTATE (-T1,V2,V3)
      CALL PRECES (T1,V3,T0,VECS)
      GO TO 50
   45 DO 48 J=1,3
   48 VECS(J) = V2(J)
C
   50 RETURN
      END
*
*   MEMBER 'VAUTL1' FOLLOWS
*
      SUBROUTINE VECTRS (RA,DEC,PMRA,PMDEC,PARLLX,RV,POS,VEL)
      save
C
C     THIS SUBROUTINE CONVERTS ANGULAR QUANTITIES TO VECTORS.
C
C          RA     = RIGHT ASCENSION IN HOURS (IN)
C          DEC    = DECLINATION IN DEGREES (IN)
C          PMRA   = PROPER MOTION IN RA IN SECONDS OF TIME PER
C                   JULIAN CENTURY (IN)
C          PMDEC  = PROPER MOTION IN DEC IN SECONDS OF ARC PER
C                   JULIAN CENTURY (IN)
C          PARLLX = PARALLAX IN SECONDS OF ARC (IN)
C          RV     = RADIAL VELOCITY IN KILOMETERS PER SECOND (IN)
C          POS    = POSITION VECTOR, EQUATORIAL RECTANGULAR COORDINATES,
C                   COMPONENTS IN AU (OUT)
C          VEL    = VELOCITY VECTOR, EQUATORIAL RECTANGULAR COORDINATES,
C                   COMPONENTS IN AU/DAY (OUT)
C
C
      DOUBLE PRECISION RA,DEC,PMRA,PMDEC,PARLLX,RV,POS,VEL,
     /     SECCON,KMAU,PARALX,DIST,R,D,CRA,SRA,CDC,SDC,PMR,PMD,RVL,
     /     DCOS,DSIN
      DIMENSION POS(3), VEL(3)
      DATA SECCON/206264.8062470964D0/,     KMAU/1.49597870D8/
C
C     IF PARALLAX IS UNKNOWN, UNDETERMINED, OR ZERO, SET IT TO 1E-7
C     SECOND OF ARC, CORRESPONDING TO A DISTANCE OF 10 MEGAPARSECS
      PARALX = PARLLX
      IF (PARALX.LE.0.0D0) PARALX = 1.0D-7
C
C     CONVERT RIGHT ASCENSION, DECLINATION, AND PARALLAX TO POSITION
C     VECTOR IN EQUATORIAL SYSTEM WITH UNITS OF AU
      DIST = SECCON / PARALX
      R = RA * 54000.0D0 / SECCON
      D = DEC * 3600.0D0 / SECCON
      CRA = DCOS(R)
      SRA = DSIN(R)
      CDC = DCOS(D)
      SDC = DSIN(D)
      POS(1) = DIST * CDC * CRA
      POS(2) = DIST * CDC * SRA
      POS(3) = DIST * SDC
C
C     CONVERT PROPER MOTION AND RADIAL VELOCITY TO ORTHOGONAL COMPONENTS
C     OF MOTION WITH UNITS OF AU/DAY
      PMR = PMRA * 15.0D0 * CDC / (PARALX * 36525.0D0)
      PMD = PMDEC / (PARALX * 36525.0D0)
      RVL = RV * 86400.0D0 / KMAU
C
C     TRANSFORM MOTION VECTOR TO EQUATORIAL SYSTEM
      VEL(1) = - PMR * SRA   - PMD * SDC * CRA   + RVL * CDC * CRA
      VEL(2) =   PMR * CRA   - PMD * SDC * SRA   + RVL * CDC * SRA
      VEL(3) =                 PMD * CDC         + RVL * SDC
C
      RETURN
      END
      SUBROUTINE ANGLES (POS,RA,DEC)
      save
C
C     THIS SUBROUTINE CONVERTS A VECTOR TO ANGULAR QUANTITIES.
C
C          POS = POSITION VECTOR, EQUATORIAL RECTANGULAR
C                COORDINATES (IN)
C          RA  = RIGHT ASCENSION IN HOURS (OUT)
C          DEC = DECLINATION IN DEGREES (OUT)
C
C
      DOUBLE PRECISION POS,RA,DEC,SECCON,XYPROJ,R,D,DSQRT,DATAN2
      DIMENSION POS(3)
      DATA SECCON/206264.8062470964D0/
      XYPROJ = DSQRT(POS(1)**2 + POS(2)**2)
      R = DATAN2(POS(2),POS(1))
      D = DATAN2(POS(3),XYPROJ)
      RA = R * SECCON / 54000.0D0
      DEC = D * SECCON / 3600.0D0
      IF (RA.LT.0.0D0) RA = RA + 24.0D0
      RETURN
      END
      SUBROUTINE PROPMO (TJD1,POS1,VEL1,TJD2,POS2)
      save
C
C     THIS SUBROUTINE APPLIES PROPER MOTION, INCLUDING FORESHORTENING
C     EFFECTS, TO A STAR'S POSITION.
C
C          TJD1 = TDB JULIAN DATE OF FIRST EPOCH (IN)
C          POS1 = POSITION VECTOR AT FIRST EPOCH (IN)
C          VEL1 = VELOCITY VECTOR AT FIRST EPOCH (IN)
C          TJD2 = TDB JULIAN DATE OF SECOND EPOCH (IN)
C          POS2 = POSITION VECTOR AT SECOND EPOCH (OUT)
C
C
      DOUBLE PRECISION TJD1,POS1,VEL1,TJD2,POS2
      DIMENSION POS1(3), VEL1(3), POS2(3)
      DO 20 J=1,3
   20 POS2(J) = POS1(J) + VEL1(J) * (TJD2 - TJD1)
      RETURN
      END
      SUBROUTINE GEOCEN (POS1,PE,POS2,TLIGHT)
      save
C
C     THIS SUBROUTINE MOVES THE ORIGIN OF COORDINATES FROM THE
C     BARYCENTER OF THE SOLAR SYSTEM TO THE CENTER OF MASS OF THE
C     EARTH, I.E., THIS SUBROUTINE CORRECTS FOR PARALLAX.
C
C          POS1   = POSITION VECTOR, REFERRED TO ORIGIN AT SOLAR SYSTEM
C                   BARYCENTER, COMPONENTS IN AU (IN)
C          PE     = POSITION VECTOR OF CENTER OF MASS OF THE EARTH,
C                   REFERRED TO ORIGIN AT SOLAR SYSTEM BARYCENTER,
C                   COMPONENTS IN AU (IN)
C          POS2   = POSITION VECTOR, REFERRED TO ORIGIN AT CENTER OF
C                   MASS OF THE EARTH, COMPONENTS IN AU (OUT)
C          TLIGHT = LIGHT TIME FROM BODY TO EARTH IN DAYS (OUT)
C
C
      DOUBLE PRECISION POS1,PE,POS2,TLIGHT,C,DSQRT
      DIMENSION POS1(3), PE(3), POS2(3)
      DATA C/173.1446333D0/
C     C = SPEED OF LIGHT IN AU/DAY = 86400 / 499.004782
      DO 20 J=1,3
   20 POS2(J) = POS1(J) - PE(J)
      TLIGHT = DSQRT(POS2(1)**2 + POS2(2)**2 + POS2(3)**2) / C
      RETURN
      END
      SUBROUTINE ABERAT (POS1,VE,TLIGHT,POS2)
      save
C
C     THIS SUBROUTINE CORRECTS POSITION VECTOR FOR ABERRATION OF LIGHT.
C     ALGORITHM INCLUDES RELATIVISTIC TERMS.  SEE MURRAY (1981)
C     MON. NOTICES ROYAL AST. SOCIETY 195, 639-648.
C
C          POS1   = POSITION VECTOR, REFERRED TO ORIGIN AT CENTER OF
C                   MASS OF THE EARTH, COMPONENTS IN AU (IN)
C          VE     = VELOCITY VECTOR OF CENTER OF MASS OF THE EARTH,
C                   REFERRED TO ORIGIN AT SOLAR SYSTEM BARYCENTER,
C                   COMPONENTS IN AU/DAY (IN)
C          TLIGHT = LIGHT TIME FROM BODY TO EARTH IN DAYS (IN)
C                   IF TLIGHT = 0.0D0, THIS SUBROUTINE WILL COMPUTE
C          POS2   = POSITION VECTOR, REFERRED TO ORIGIN AT CENTER OF
C                   MASS OF THE EARTH, CORRECTED FOR ABERRATION,
C                   COMPONENTS IN AU (OUT)
C
C
      DOUBLE PRECISION POS1,VE,TLIGHT,POS2,C,TL,P1MAG,VEMAG,
     /     BETA,DOT,COSD,GAMMAI,P,Q,R,DSQRT
      DIMENSION POS1(3), VE(3), POS2(3)
      DATA C/173.1446333D0/
C     C = SPEED OF LIGHT IN AU/DAY = 86400 / 499.004782
C
      TL = TLIGHT
      P1MAG = TL * C
      IF (TL.NE.0.0D0) GO TO 20
      P1MAG = DSQRT(POS1(1)**2 + POS1(2)**2 + POS1(3)**2)
      TL = P1MAG / C
   20 VEMAG = DSQRT(VE(1)**2 + VE(2)**2 + VE(3)**2)
      BETA = VEMAG / C
      DOT = POS1(1)*VE(1) + POS1(2)*VE(2) + POS1(3)*VE(3)
      COSD = DOT / (P1MAG * VEMAG)
      GAMMAI = DSQRT(1.0D0 - BETA**2)
      P = BETA * COSD
      Q = (1.0D0 + P / (1.0D0 + GAMMAI)) * TL
      R = 1.0D0 + P
C
      DO 30 J=1,3
   30 POS2(J) = (GAMMAI * POS1(J) + Q * VE(J)) / R
      RETURN
      END
      SUBROUTINE PRECES (TJD1,POS1,TJD2,POS2)
      save
C
C     THIS SUBROUTINE PRECESSES EQUATORIAL RECTANGULAR COORDINATES FROM
C     ONE EPOCH TO ANOTHER.  THE COORDINATES ARE REFERRED TO THE MEAN
C     EQUATOR AND EQUINOX OF THE TWO RESPECTIVE EPOCHS.  SEE PAGES 30-34
C     OF THE EXPLANATORY SUPPLEMENT TO THE AE, LIESKE, ET AL. (1977)
C     ASTRONOMY AND ASTROPHYSICS 58, 1-16, AND LIESKE (1979) ASTRONOMY
C     AND ASTROPHYSICS 73, 282-284.
C
C          TJD1 = TDB JULIAN DATE OF FIRST EPOCH (IN)
C          POS1 = POSITION VECTOR, GEOCENTRIC EQUATORIAL RECTANGULAR
C                 COORDINATES, REFERRED TO MEAN EQUATOR AND EQUINOX OF
C                 FIRST EPOCH (IN)
C          TJD2 = TDB JULIAN DATE OF SECOND EPOCH (IN)
C          POS2 = POSITION VECTOR, GEOCENTRIC EQUATORIAL RECTANGULAR
C                 COORDINATES, REFERRED TO MEAN EQUATOR AND EQUINOX OF
C                 SECOND EPOCH (OUT)
C
C
      DOUBLE PRECISION TJD1,TJD2,T0,T,T02,T2,T3,POS1,POS2,SECCON,
     /     ZETA0,ZEE,THETA,CZETA0,SZETA0,CZEE,SZEE,CTHETA,STHETA,
     /     XX,YX,ZX,XY,YY,ZY,XZ,YZ,ZZ,T1LAST,T2LAST,DABS,DCOS,DSIN
      DIMENSION POS1(3), POS2(3)
      DATA SECCON/206264.8062470964D0/
      DATA T1LAST,T2LAST/0.0D0,0.0D0/
C
      IF (DABS(TJD1-T1LAST).LT.1.0D-6.AND.DABS(TJD2-T2LAST).LT.1.0D-6)
     /     GO TO 20
      IF (DABS(TJD1-T2LAST).LT.1.0D-6.AND.DABS(TJD2-T1LAST).LT.1.0D-6)
     /     GO TO 30
C
C     T0 AND T BELOW CORRESPOND TO LIESKE'S BIG T AND LITTLE T
      T0 = (TJD1 - 2451545.0D0) / 36525.0D0
      T = (TJD2 - TJD1) / 36525.0D0
      T02 = T0 * T0
      T2 = T * T
      T3 = T2 * T
C     ZETA0, ZEE, AND THETA BELOW CORRESPOND TO LIESKE'S ZETA-SUB-A,
C     Z-SUB-A, AND THETA-SUB-A
      ZETA0 = (2306.2181D0 + 1.39656D0*T0 - 0.000139D0*T02) * T
     /      + (0.30188D0 - 0.000344D0*T0) * T2
     /      +  0.017998D0 * T3
      ZEE   = (2306.2181D0 + 1.39656D0*T0 - 0.000139D0*T02) * T
     /      + (1.09468D0 + 0.000066D0*T0) * T2
     /      +  0.018203D0 * T3
      THETA = (2004.3109D0 - 0.85330D0*T0 - 0.000217D0*T02) * T
     /      + (-0.42665D0 - 0.000217D0*T0) * T2
     /      -  0.041833D0 * T3
      ZETA0 = ZETA0 / SECCON
      ZEE = ZEE / SECCON
      THETA = THETA / SECCON
      CZETA0 = DCOS(ZETA0)
      SZETA0 = DSIN(ZETA0)
      CZEE = DCOS(ZEE)
      SZEE = DSIN(ZEE)
      CTHETA = DCOS(THETA)
      STHETA = DSIN(THETA)
C
C     PRECESSION ROTATION MATRIX FOLLOWS
      XX = CZETA0*CTHETA*CZEE - SZETA0*SZEE
      YX = -SZETA0*CTHETA*CZEE - CZETA0*SZEE
      ZX = -STHETA*CZEE
      XY = CZETA0*CTHETA*SZEE + SZETA0*CZEE
      YY = -SZETA0*CTHETA*SZEE + CZETA0*CZEE
      ZY = -STHETA*SZEE
      XZ = CZETA0*STHETA
      YZ = -SZETA0*STHETA
      ZZ = CTHETA
      T1LAST = TJD1
      T2LAST = TJD2
C
C     PERFORM ROTATION
   20 POS2(1) = XX*POS1(1) + YX*POS1(2) + ZX*POS1(3)
      POS2(2) = XY*POS1(1) + YY*POS1(2) + ZY*POS1(3)
      POS2(3) = XZ*POS1(1) + YZ*POS1(2) + ZZ*POS1(3)
      GO TO 50
C
C     PERFORM INVERSE ROTATION
   30 POS2(1) = XX*POS1(1) + XY*POS1(2) + XZ*POS1(3)
      POS2(2) = YX*POS1(1) + YY*POS1(2) + YZ*POS1(3)
      POS2(3) = ZX*POS1(1) + ZY*POS1(2) + ZZ*POS1(3)
C
   50 RETURN
      END
      SUBROUTINE NUTATE (TJD,POS1,POS2)
      save
C
C     THIS SUBROUTINE NUTATES EQUATORIAL RECTANGULAR COORDINATES FROM
C     MEAN EQUATOR AND EQUINOX OF EPOCH TO TRUE EQUATOR AND EQUINOX OF
C     EPOCH.  SEE PAGES 41-45 OF THE EXPLANATORY SUPPLEMENT TO THE AE.
C
C          TJD    = TDB JULIAN DATE OF EPOCH (IN)
C          POS1   = POSITION VECTOR, GEOCENTRIC EQUATORIAL RECTANGULAR
C                   COORDINATES, REFERRED TO MEAN EQUATOR AND EQUINOX
C                   OF EPOCH (IN)
C          POS2   = POSITION VECTOR, GEOCENTRIC EQUATORIAL RECTANGULAR
C                   COORDINATES, REFERRED TO TRUE EQUATOR AND EQUINOX
C                   OF EPOCH (OUT)
C
C          NOTE --  IF TJD IS NEGATIVE, INVERSE NUTATION (TRUE TO MEAN)
C                   IS APPLIED
C
C
      DOUBLE PRECISION TJD,POS1,POS2,TLAST,TJD1,SECCON,OBLM,OBLT,EQEQ,
     /     DPSI,DEPS,COBM,SOBM,COBT,SOBT,CPSI,SPSI,
     /     XX,YX,ZX,XY,YY,ZY,XZ,YZ,ZZ,DABS,DCOS,DSIN
      DIMENSION POS1(3), POS2(3)
      DATA SECCON/206264.8062470964D0/
      DATA TLAST/0.0D0/
C
      TJD1 = DABS(TJD)
      IF (DABS(TJD1-TLAST).LT.1.0D-6) GO TO 10
C
      CALL ETILT (TJD1,OBLM,OBLT,EQEQ,DPSI,DEPS)
      OBLM = OBLM * 3600.0D0 / SECCON
      OBLT = OBLT * 3600.0D0 / SECCON
      DPSI = DPSI / SECCON
      DEPS = DEPS / SECCON
      COBM = DCOS(OBLM)
      SOBM = DSIN(OBLM)
      COBT = DCOS(OBLT)
      SOBT = DSIN(OBLT)
      CPSI = DCOS(DPSI)
      SPSI = DSIN(DPSI)
C
C     NUTATION ROTATION MATRIX FOLLOWS
      XX = CPSI
      YX = -SPSI*COBM
      ZX = -SPSI*SOBM
      XY = SPSI*COBT
      YY = CPSI*COBM*COBT + SOBM*SOBT
      ZY = CPSI*SOBM*COBT - COBM*SOBT
      XZ = SPSI*SOBT
      YZ = CPSI*COBM*SOBT - SOBM*COBT
      ZZ = CPSI*SOBM*SOBT + COBM*COBT
      TLAST = TJD1
   10 IF (TJD.LT.0.0D0) GO TO 30
C
C     PERFORM ROTATION
   20 POS2(1) = XX*POS1(1) + YX*POS1(2) + ZX*POS1(3)
      POS2(2) = XY*POS1(1) + YY*POS1(2) + ZY*POS1(3)
      POS2(3) = XZ*POS1(1) + YZ*POS1(2) + ZZ*POS1(3)
      GO TO 50
C
C     PERFORM INVERSE ROTATION
   30 POS2(1) = XX*POS1(1) + XY*POS1(2) + XZ*POS1(3)
      POS2(2) = YX*POS1(1) + YY*POS1(2) + YZ*POS1(3)
      POS2(3) = ZX*POS1(1) + ZY*POS1(2) + ZZ*POS1(3)
C
   50 RETURN
      END
      SUBROUTINE SPIN (ST,POS1,POS2)
      save
C
C     THIS SUBROUTINE TRANSFORMS GEOCENTRIC RECTANGULAR COORDINATES
C     FROM ROTATING SYSTEM BASED ON ROTATIONAL EQUATOR AND ORTHOGONAL
C     REFERENCE MERIDIAN TO NON-ROTATING SYSTEM BASED ON TRUE EQUATOR
C     AND EQUINOX OF DATE.
C
C          ST     = LOCAL APPARENT SIDEREAL TIME AT REFERENCE MERIDIAN
C                   IN HOURS (IN)
C          POS1   = VECTOR IN GEOCENTRIC RECTANGULAR
C                   ROTATING SYSTEM, REFERRED TO ROTATIONAL EQUATOR
C                   AND ORTHOGONAL REFERENCE MERIDIAN (IN)
C          POS2   = VECTOR IN GEOCENTRIC RECTANGULAR
C                   NON-ROTATING SYSTEM, REFERRED TO TRUE EQUATOR
C                   AND EQUINOX OF DATE (OUT)
C
C
      DOUBLE PRECISION ST,POS1,POS2,SECCON,TLAST,STR,COSST,SINST,
     /     XX,YX,ZX,XY,YY,ZY,XZ,YZ,ZZ,DABS,DCOS,DSIN
      DIMENSION POS1(3), POS2(3)
C
      DATA SECCON/206264.8062470964D0/
      DATA TLAST/-999.0D0/
C
      IF (DABS(ST-TLAST).LT.1.0D-12) GO TO 10
C
      STR   = ST * 15.0D0 * 3600.0D0 / SECCON
      COSST = DCOS(STR)
      SINST = DSIN(STR)
C
C     SIDEREAL TIME ROTATION MATRIX FOLLOWS
      XX =  COSST
      YX = -SINST
      ZX =  0.0D0
      XY =  SINST
      YY =  COSST
      ZY =  0.0D0
      XZ =  0.0D0
      YZ =  0.0D0
      ZZ =  1.0D0
      TLAST = ST
   10 CONTINUE
C
C     PERFORM ROTATION
   20 POS2(1) = XX*POS1(1) + YX*POS1(2) + ZX*POS1(3)
      POS2(2) = XY*POS1(1) + YY*POS1(2) + ZY*POS1(3)
      POS2(3) = XZ*POS1(1) + YZ*POS1(2) + ZZ*POS1(3)
C
   50 RETURN
      END
      SUBROUTINE WOBBLE (X,Y,POS1,POS2)
      save
C
C     THIS SUBROUTINE CORRECTS EARTH-FIXED GEOCENTRIC RECTANGULAR
C     COORDINATES FOR POLAR MOTION.  IT TRANSFORMS A VECTOR FROM
C     EARTH-FIXED GEOGRAPHIC SYSTEM TO ROTATING SYSTEM BASED ON
C     ROTATIONAL EQUATOR AND ORTHOGONAL GREENWICH MERIDIAN THROUGH
C     AXIS OF ROTATION.
C
C          X      = CONVENTIONALLY-DEFINED X COORDINATE OF ROTATIONAL
C                   POLE WITH RESPECT TO CIO, IN ARCSECONDS (IN)
C          Y      = CONVENTIONALLY-DEFINED Y COORDINATE OF ROTATIONAL
C                   POLE WITH RESPECT TO CIO, IN ARCSECONDS (IN)
C          POS1   = VECTOR IN GEOCENTRIC RECTANGULAR
C                   EARTH-FIXED SYSTEM, REFERRED TO GEOGRAPHIC
C                   EQUATOR AND GREENWICH MERIDIAN (IN)
C          POS2   = VECTOR IN GEOCENTRIC RECTANGULAR
C                   ROTATING SYSTEM, REFERRED TO ROTATIONAL EQUATOR
C                   AND ORTHOGONAL GREENWICH MERIDIAN (OUT)
C
C
      DOUBLE PRECISION X,Y,POS1,POS2,SECCON,XPOLE,YPOLE,
     /     XX,YX,ZX,XY,YY,ZY,XZ,YZ,ZZ
      DIMENSION POS1(3), POS2(3)
C
      DATA SECCON/206264.8062470964D0/
C
      XPOLE = X / SECCON
      YPOLE = Y / SECCON
C
C     WOBBLE ROTATION MATRIX FOLLOWS
      XX =  1.0D0
      YX =  0.0D0
      ZX = -XPOLE
      XY =  0.0D0
      YY =  1.0D0
      ZY =  YPOLE
      XZ =  XPOLE
      YZ = -YPOLE
      ZZ =  1.0D0
   10 CONTINUE
C
C     PERFORM ROTATION
   20 POS2(1) = XX*POS1(1) + YX*POS1(2) + ZX*POS1(3)
      POS2(2) = XY*POS1(1) + YY*POS1(2) + ZY*POS1(3)
      POS2(3) = XZ*POS1(1) + YZ*POS1(2) + ZZ*POS1(3)
C
   50 RETURN
      END
      SUBROUTINE TERRA (GLON,GLAT,HT,ST,POS,VEL)
      save
C
C     THIS SUBROUTINE COMPUTES THE POSITION AND VELOCITY VECTORS OF
C     A TERRESTRIAL OBSERVER WITH RESPECT TO THE CENTER OF THE EARTH.
C
C          GLON   = LONGITUDE OF OBSERVER WITH RESPECT TO REFERENCE
C                   MERIDIAN (EAST +) IN DEGREES (IN)
C          GLAT   = GEODETIC LATITUDE (NORTH +) OF OBSERVER
C                   IN DEGREES (IN)
C          HT     = HEIGHT OF OBSERVER IN METERS (IN)
C          ST     = LOCAL APPARENT SIDEREAL TIME AT REFERENCE MERIDIAN
C                   IN HOURS (IN)
C          POS    = POSITION VECTOR OF OBSERVER WITH RESPECT TO CENTER
C                   OF EARTH, EQUATORIAL RECTANGULAR COORDINATES,
C                   REFERRED TO TRUE EQUATOR AND EQUINOX OF DATE,
C                   COMPONENTS IN AU (OUT)
C          VEL    = VELOCITY VECTOR OF OBSERVER WITH RESPECT TO CENTER
C                   OF EARTH, EQUATORIAL RECTANGULAR COORDINATES,
C                   REFERRED TO TRUE EQUATOR AND EQUINOX OF DATE,
C                   COMPONENTS IN AU/DAY (OUT)
C
C          NOTE --  IF REFERENCE MERIDIAN IS GREENWICH AND ST=0, POS
C                   IS EFFECTIVELY REFERRED TO EQUATOR AND GREENWICH
C
C
      DOUBLE PRECISION GLON,GLAT,HT,ST,POS,VEL,SECCON,ERAD,F,OMEGA,
     /     KMAU,DF2,PHI,SINPHI,COSPHI,C,S,ACH,ASH,STLOCL,SINST,COSST,
     /     DSQRT,DCOS,DSIN
      DIMENSION POS(3), VEL(3)
C
      DATA SECCON/206264.8062470964D0/
C
      DATA ERAD / 6378.140D0 /,   F / 0.00335281D0 /
C     ERAD = RADIUS OF EARTH IN KM, F = EARTH ELLIPSOID FLATTENING
      DATA OMEGA / 7.2921151467D-5 /
C     OMEGA = ROTATIONAL ANGULAR VELOCITY OF EARTH IN RADIANS/SEC
      DATA KMAU / 1.49597870D8 /
C     KMAU = KILOMETERS PER ASTRONOMICAL UNIT
C
C     COMPUTE PARAMETERS RELATING TO GEODETIC TO GEOCENTRIC CONVERSION
      DF2 = (1.0D0 - F)**2
      PHI = GLAT * 3600.0D0 / SECCON
      SINPHI = DSIN(PHI)
      COSPHI = DCOS(PHI)
      C = 1.0D0 / DSQRT ( COSPHI**2 + DF2 * SINPHI**2 )
      S = DF2 * C
      ACH = ERAD * C + HT/1000.0D0
      ASH = ERAD * S + HT/1000.0D0
C
C     COMPUTE LOCAL SIDEREAL TIME FACTORS
      STLOCL = (ST * 54000.0D0 + GLON * 3600.0D0) / SECCON
      SINST = DSIN(STLOCL)
      COSST = DCOS(STLOCL)
C
C     COMPUTE POSITION VECTOR COMPONENTS IN KM
      POS(1) = ACH * COSPHI * COSST
      POS(2) = ACH * COSPHI * SINST
      POS(3) = ASH * SINPHI
C
C     COMPUTE VELOCITY VECTOR COMPONENTS IN KM/SEC
      VEL(1) = -OMEGA * ACH * COSPHI * SINST
      VEL(2) =  OMEGA * ACH * COSPHI * COSST
      VEL(3) =  0.0D0
C
C     CONVERT POSITION AND VELOCITY COMPONENTS TO AU AND AU/DAY
      DO 20 J=1,3
      POS(J) = POS(J) / KMAU
      VEL(J) = VEL(J) / KMAU * 86400.0D0
   20 CONTINUE
C
      RETURN
      END
*
*   MEMBER 'VABAS1' FOLLOWS
*
      SUBROUTINE TIMES (TDBJD,TDTJD,SECDIF)
      save
C
C     THIS SUBROUTINE COMPUTES THE TERRESTRIAL DYNAMICAL
C     TIME (TDT) JULIAN DATE CORRESPONDING TO A BARYCENTRIC
C     DYNAMICAL TIME (TDB) JULIAN DATE.
C     EXPRESSIONS USED IN THIS VERSION ARE APPROXIMATIONS
C     RESULTING IN ACCURACIES OF ABOUT 20 MICROSECONDS.
C
C          TDBJD  = TDB JULIAN DATE (IN)
C          TDTJD  = TDT JULIAN DATE (OUT)
C          SECDIF = DIFFERENCE TDBJD-TDTJD, IN SECONDS (OUT)
C
C
      DOUBLE PRECISION TDBJD,TDTJD,SECDIF,T,M,LLJ,E,DSIN
C
      T = (TDBJD - 2433282.5D0) / 36525.0D0
      M = 6.248291D0 + 628.3019415D0 * T
      LLJ = 5.652593D0 + 575.3380832D0 * T
      E = M + 0.01671D0 * DSIN(M)
      SECDIF = 1.658D-3 * DSIN(E) + 20.73D-6 * DSIN(LLJ)
      TDTJD = TDBJD - SECDIF / 86400.0D0
C
      RETURN
      END
      SUBROUTINE ETILT (TJD,OBLM,OBLT,EQEQ,DPSI,DEPS)
      save
C
C     THIS SUBROUTINE COMPUTES QUANTITIES RELATED TO THE ORIENTATION
C     OF THE EARTH'S ROTATION AXIS AT JULIAN DATE TJD.
C
C          TJD    = TDB JULIAN DATE FOR ORIENTATION PARAMETERS (IN)
C          OBLM   = MEAN OBLIQUITY OF THE ECLIPTIC IN DEGREES AT
C                   DATE TJD (OUT)
C          OBLT   = TRUE OBLIQUITY OF THE ECLIPTIC IN DEGREES AT
C                   DATE TJD (OUT)
C          EQEQ   = EQUATION OF THE EQUINOXES IN SECONDS OF TIME AT
C                   DATE TJD (OUT)
C          DPSI   = NUTATION IN LONGITUDE IN SECONDS OF ARC AT
C                   DATE TJD (OUT)
C          DEPS   = NUTATION IN OBLIQUITY IN SECONDS OF ARC AT
C                   DATE TJD (OUT)
C
C
      DOUBLE PRECISION TJD,T0,T,T2,T3,TLAST,OBLM,OBLT,EQEQ,DPSI,DEPS,
     /     SECCON,OBM,OBT,EE,PSI,EPS,DABS,DCOS
C
      DATA T0/2451545.0D0/
C     T0 = TDB JULIAN DATE OF EPOCH J2000.0
      DATA SECCON/206264.8062470964D0/
      DATA TLAST/0.0D0/
C
      IF (DABS(TJD-TLAST).LT.1.0D-6) GO TO 20
      T = (TJD - T0) / 36525.0D0
      T2 = T * T
      T3 = T2 * T
C
C     OBTAIN NUTATION PARAMETERS IN SECONDS OF ARC
      CALL NOD (T,PSI,EPS)
C
C     COMPUTE MEAN OBLIQUITY OF THE ECLIPTIC IN SECONDS OF ARC
      OBM = 84381.4480D0 - 46.8150D0*T - 0.00059D0*T2
     /    + 0.001813D0*T3
C
C     COMPUTE TRUE OBLIQUITY OF THE ECLIPTIC IN SECONDS OF ARC
      OBT = OBM + EPS
C
C     COMPUTE EQUATION OF THE EQUINOXES IN SECONDS OF TIME
      EE = PSI / 15.0D0 * DCOS (OBT/SECCON)
C
C     CONVERT OBLIQUITY VALUES TO DEGREES
      OBM = OBM / 3600.0D0
      OBT = OBT / 3600.0D0
      TLAST = TJD
C
   20 OBLM = OBM
      OBLT = OBT
      EQEQ = EE
      DPSI = PSI
      DEPS = EPS
C
      RETURN
      END
*
*   MEMBER 'VASUN1' FOLLOWS
*
      SUBROUTINE SUNFLD (POS1,PE,POS2)
      save
C
C     THIS SUBROUTINE CORRECTS POSITION VECTOR FOR THE DEFLECTION
C     OF LIGHT IN THE GRAVITATIONAL FIELD OF THE SUN.  SEE MISNER,
C     THORNE, AND WHEELER (1973), GRAVITATION, PP. 184-185.  THIS
C     SUBROUTINE VALID FOR BODIES WITHIN THE SOLAR SYSTEM AS WELL AS
C     FOR STARS.
C
C          POS1 = POSITION VECTOR, REFERRED TO ORIGIN AT CENTER OF MASS
C                 OF THE EARTH, COMPONENTS IN AU (IN)
C          PE   = POSITION VECTOR OF CENTER OF MASS OF THE EARTH,
C                 REFERRED TO ORIGIN AT CENTER OF MASS OF
C                 THE SUN, COMPONENTS IN AU (IN)
C          POS2 = POSITION VECTOR, REFERRED TO ORIGIN AT CENTER OF MASS
C                 OF THE EARTH, CORRECTED FOR GRAVITATIONAL DEFLEC-
C                 TION, COMPONENTS IN AU (OUT)
C
C
      DOUBLE PRECISION POS1,PE,POS2,P1HAT,PEHAT,P1MAG,PEMAG,MAU,GS,C,F,
     /     COSD,SIND,B,BM,PQMAG,ZFINL,ZINIT,XIFINL,XIINIT,
     /     DELPHI,DELPHP,DELP,DABS,DSQRT
      DIMENSION POS1(3), PE(3), POS2(3), P1HAT(3), PEHAT(3)
      DATA MAU/1.49597870D11/
C     MAU = NUMBER OF METERS PER AU
      DATA GS/1.32712438D20/
C     GS = HELIOCENTRIC GRAVITATIONAL CONSTANT
      DATA C/299792458.0D0/
C     C = SPEED OF LIGHT
      F = 0.0D0
C
C     COMPUTE VECTOR MAGNITUDES AND UNIT VECTORS
      P1MAG = DSQRT (POS1(1)**2 + POS1(2)**2 + POS1(3)**2)
      PEMAG = DSQRT (  PE(1)**2 +   PE(2)**2 +   PE(3)**2)
      DO 20 J=1,3
      P1HAT(J) = POS1(J) / P1MAG
   20 PEHAT(J) =   PE(J) / PEMAG
C
C     COMPUTE GEOMETRICAL QUANTITIES
C     COSD AND SIND ARE COSINE AND SINE OF D, THE ANGULAR SEPARATION
C     OF THE BODY FROM THE SUN AS VIEWED FROM THE EARTH
      COSD = - PEHAT(1)*P1HAT(1) - PEHAT(2)*P1HAT(2) - PEHAT(3)*P1HAT(3)
      IF (DABS(COSD).GT.0.9999999999D0) GO TO 40
      SIND = DSQRT (1.0D0 - COSD**2)
C     B IS THE IMPACT PARAMETER FOR THE RAY
      B = PEMAG * SIND
      BM = B * MAU
C     PQMAG IS THE DISTANCE OF THE BODY FROM THE SUN
      PQMAG = DSQRT (P1MAG**2 + PEMAG**2 - 2.0D0 * P1MAG * PEMAG * COSD)
C
C     COMPUTE DELPHI, THE ANGLE OF DEFLECTION OF THE RAY
      ZFINL = PEMAG * COSD
      ZINIT = -P1MAG + ZFINL
      XIFINL = ZFINL / B
      XIINIT = ZINIT / B
      DELPHI = 2.0D0*GS/(BM*C*C) * (XIFINL / DSQRT (1.0D0 + XIFINL**2)
     /                            - XIINIT / DSQRT (1.0D0 + XIINIT**2))
C
C     COMPUTE DELPHP, THE CHANGE IN ANGLE AS SEEN AT THE EARTH
      DELPHP = DELPHI / (1.0D0 + (PEMAG / PQMAG))
C
C     FIX UP POSITION VECTOR
C     POS2 IS POS1 ROTATED THROUGH ANGLE DELPHP IN PLANE DEFINED
C     BY POS1 AND PE
      F = DELPHP * P1MAG / SIND
   40 DO 50 J=1,3
      DELP = F * (COSD * P1HAT(J) + PEHAT(J))
   50 POS2(J) = POS1(J) + DELP
C
      RETURN
      END
*
*   MEMBER 'VANUT1' FOLLOWS
*
      SUBROUTINE NOD (T,DPSI,DEPS)
      save
C
C     THIS SUBROUTINE EVALUATES THE NUTATION SERIES AND RETURNS THE
C     VALUES FOR NUTATION IN LONGITUDE AND NUTATION IN OBLIQUITY.
C     WAHR NUTATION SERIES FOR AXIS B FOR GILBERT & DZIEWONSKI EARTH
C     MODEL 1066A.
C     1980 IAU THEORY OF NUTATION.
C
C          T    = TDB TIME IN JULIAN CENTURIES SINCE J2000.0 (IN)
C          DPSI = NUTATION IN LONGITUDE IN SECONDS OF ARC (OUT)
C          DEPS = NUTATION IN OBLIQUITY IN SECONDS OF ARC (OUT)
C
C
      DOUBLE PRECISION T,DPSI,DEPS,SECCON,L,LP,F,D,OM,ARG,DMOD,DBLE,
     /     DSIN,DCOS
      DIMENSION X(9,106),X1(90),X2(90),X3(90),X4(90),X5(90),X6(90),
     /     X7(90),X8(90),X9(90),XA(90),XB(54)
      EQUIVALENCE(X(1,  1),X1(1))
      EQUIVALENCE(X(1, 11),X2(1))
      EQUIVALENCE(X(1, 21),X3(1))
      EQUIVALENCE(X(1, 31),X4(1))
      EQUIVALENCE(X(1, 41),X5(1))
      EQUIVALENCE(X(1, 51),X6(1))
      EQUIVALENCE(X(1, 61),X7(1))
      EQUIVALENCE(X(1, 71),X8(1))
      EQUIVALENCE(X(1, 81),X9(1))
      EQUIVALENCE(X(1, 91),XA(1))
      EQUIVALENCE(X(1,101),XB(1))
      DATA SECCON/206264.8062470964D0/
C
C***********************************************************************
C
C
C     TABLE OF MULTIPLES OF ARGUMENTS AND COEFFICIENTS
C
C                   MULTIPLE OF            LONGITUDE        OBLIQUITY
C              L    L'   F    D  OMEGA   COEFF. OF SIN    COEFF. OF COS
      DATA X1/ 0.,  0.,  0.,  0.,  1., -171996., -174.2,  92025.,  8.9,
     /         0.,  0.,  2., -2.,  2.,  -13187.,   -1.6,   5736., -3.1,
     /         0.,  0.,  2.,  0.,  2.,   -2274.,   -0.2,    977., -0.5,
     /         0.,  0.,  0.,  0.,  2.,    2062.,    0.2,   -895.,  0.5,
     /         0.,  1.,  0.,  0.,  0.,    1426.,   -3.4,     54., -0.1,
     /         1.,  0.,  0.,  0.,  0.,     712.,    0.1,     -7.,  0.0,
     /         0.,  1.,  2., -2.,  2.,    -517.,    1.2,    224., -0.6,
     /         0.,  0.,  2.,  0.,  1.,    -386.,   -0.4,    200.,  0.0,
     /         1.,  0.,  2.,  0.,  2.,    -301.,    0.0,    129., -0.1,
     /         0., -1.,  2., -2.,  2.,     217.,   -0.5,    -95.,  0.3/
      DATA X2/ 1.,  0.,  0., -2.,  0.,    -158.,    0.0,     -1.,  0.0,
     /         0.,  0.,  2., -2.,  1.,     129.,    0.1,    -70.,  0.0,
     /        -1.,  0.,  2.,  0.,  2.,     123.,    0.0,    -53.,  0.0,
     /         1.,  0.,  0.,  0.,  1.,      63.,    0.1,    -33.,  0.0,
     /         0.,  0.,  0.,  2.,  0.,      63.,    0.0,     -2.,  0.0,
     /        -1.,  0.,  2.,  2.,  2.,     -59.,    0.0,     26.,  0.0,
     /        -1.,  0.,  0.,  0.,  1.,     -58.,   -0.1,     32.,  0.0,
     /         1.,  0.,  2.,  0.,  1.,     -51.,    0.0,     27.,  0.0,
     /         2.,  0.,  0., -2.,  0.,      48.,    0.0,      1.,  0.0,
     /        -2.,  0.,  2.,  0.,  1.,      46.,    0.0,    -24.,  0.0/
      DATA X3/ 0.,  0.,  2.,  2.,  2.,     -38.,    0.0,     16.,  0.0,
     /         2.,  0.,  2.,  0.,  2.,     -31.,    0.0,     13.,  0.0,
     /         2.,  0.,  0.,  0.,  0.,      29.,    0.0,     -1.,  0.0,
     /         1.,  0.,  2., -2.,  2.,      29.,    0.0,    -12.,  0.0,
     /         0.,  0.,  2.,  0.,  0.,      26.,    0.0,     -1.,  0.0,
     /         0.,  0.,  2., -2.,  0.,     -22.,    0.0,      0.,  0.0,
     /        -1.,  0.,  2.,  0.,  1.,      21.,    0.0,    -10.,  0.0,
     /         0.,  2.,  0.,  0.,  0.,      17.,   -0.1,      0.,  0.0,
     /         0.,  2.,  2., -2.,  2.,     -16.,    0.1,      7.,  0.0,
     /        -1.,  0.,  0.,  2.,  1.,      16.,    0.0,     -8.,  0.0/
      DATA X4/ 0.,  1.,  0.,  0.,  1.,     -15.,    0.0,      9.,  0.0,
     /         1.,  0.,  0., -2.,  1.,     -13.,    0.0,      7.,  0.0,
     /         0., -1.,  0.,  0.,  1.,     -12.,    0.0,      6.,  0.0,
     /         2.,  0., -2.,  0.,  0.,      11.,    0.0,      0.,  0.0,
     /        -1.,  0.,  2.,  2.,  1.,     -10.,    0.0,      5.,  0.0,
     /         1.,  0.,  2.,  2.,  2.,      -8.,    0.0,      3.,  0.0,
     /         0., -1.,  2.,  0.,  2.,      -7.,    0.0,      3.,  0.0,
     /         0.,  0.,  2.,  2.,  1.,      -7.,    0.0,      3.,  0.0,
     /         1.,  1.,  0., -2.,  0.,      -7.,    0.0,      0.,  0.0,
     /         0.,  1.,  2.,  0.,  2.,       7.,    0.0,     -3.,  0.0/
      DATA X5/-2.,  0.,  0.,  2.,  1.,      -6.,    0.0,      3.,  0.0,
     /         0.,  0.,  0.,  2.,  1.,      -6.,    0.0,      3.,  0.0,
     /         2.,  0.,  2., -2.,  2.,       6.,    0.0,     -3.,  0.0,
     /         1.,  0.,  0.,  2.,  0.,       6.,    0.0,      0.,  0.0,
     /         1.,  0.,  2., -2.,  1.,       6.,    0.0,     -3.,  0.0,
     /         0.,  0.,  0., -2.,  1.,      -5.,    0.0,      3.,  0.0,
     /         0., -1.,  2., -2.,  1.,      -5.,    0.0,      3.,  0.0,
     /         2.,  0.,  2.,  0.,  1.,      -5.,    0.0,      3.,  0.0,
     /         1., -1.,  0.,  0.,  0.,       5.,    0.0,      0.,  0.0,
     /         1.,  0.,  0., -1.,  0.,      -4.,    0.0,      0.,  0.0/
      DATA X6/ 0.,  0.,  0.,  1.,  0.,      -4.,    0.0,      0.,  0.0,
     /         0.,  1.,  0., -2.,  0.,      -4.,    0.0,      0.,  0.0,
     /         1.,  0., -2.,  0.,  0.,       4.,    0.0,      0.,  0.0,
     /         2.,  0.,  0., -2.,  1.,       4.,    0.0,     -2.,  0.0,
     /         0.,  1.,  2., -2.,  1.,       4.,    0.0,     -2.,  0.0,
     /         1.,  1.,  0.,  0.,  0.,      -3.,    0.0,      0.,  0.0,
     /         1., -1.,  0., -1.,  0.,      -3.,    0.0,      0.,  0.0,
     /        -1., -1.,  2.,  2.,  2.,      -3.,    0.0,      1.,  0.0,
     /         0., -1.,  2.,  2.,  2.,      -3.,    0.0,      1.,  0.0,
     /         1., -1.,  2.,  0.,  2.,      -3.,    0.0,      1.,  0.0/
      DATA X7/ 3.,  0.,  2.,  0.,  2.,      -3.,    0.0,      1.,  0.0,
     /        -2.,  0.,  2.,  0.,  2.,      -3.,    0.0,      1.,  0.0,
     /         1.,  0.,  2.,  0.,  0.,       3.,    0.0,      0.,  0.0,
     /        -1.,  0.,  2.,  4.,  2.,      -2.,    0.0,      1.,  0.0,
     /         1.,  0.,  0.,  0.,  2.,      -2.,    0.0,      1.,  0.0,
     /        -1.,  0.,  2., -2.,  1.,      -2.,    0.0,      1.,  0.0,
     /         0., -2.,  2., -2.,  1.,      -2.,    0.0,      1.,  0.0,
     /        -2.,  0.,  0.,  0.,  1.,      -2.,    0.0,      1.,  0.0,
     /         2.,  0.,  0.,  0.,  1.,       2.,    0.0,     -1.,  0.0,
     /         3.,  0.,  0.,  0.,  0.,       2.,    0.0,      0.,  0.0/
      DATA X8/ 1.,  1.,  2.,  0.,  2.,       2.,    0.0,     -1.,  0.0,
     /         0.,  0.,  2.,  1.,  2.,       2.,    0.0,     -1.,  0.0,
     /         1.,  0.,  0.,  2.,  1.,      -1.,    0.0,      0.,  0.0,
     /         1.,  0.,  2.,  2.,  1.,      -1.,    0.0,      1.,  0.0,
     /         1.,  1.,  0., -2.,  1.,      -1.,    0.0,      0.,  0.0,
     /         0.,  1.,  0.,  2.,  0.,      -1.,    0.0,      0.,  0.0,
     /         0.,  1.,  2., -2.,  0.,      -1.,    0.0,      0.,  0.0,
     /         0.,  1., -2.,  2.,  0.,      -1.,    0.0,      0.,  0.0,
     /         1.,  0., -2.,  2.,  0.,      -1.,    0.0,      0.,  0.0,
     /         1.,  0., -2., -2.,  0.,      -1.,    0.0,      0.,  0.0/
      DATA X9/ 1.,  0.,  2., -2.,  0.,      -1.,    0.0,      0.,  0.0,
     /         1.,  0.,  0., -4.,  0.,      -1.,    0.0,      0.,  0.0,
     /         2.,  0.,  0., -4.,  0.,      -1.,    0.0,      0.,  0.0,
     /         0.,  0.,  2.,  4.,  2.,      -1.,    0.0,      0.,  0.0,
     /         0.,  0.,  2., -1.,  2.,      -1.,    0.0,      0.,  0.0,
     /        -2.,  0.,  2.,  4.,  2.,      -1.,    0.0,      1.,  0.0,
     /         2.,  0.,  2.,  2.,  2.,      -1.,    0.0,      0.,  0.0,
     /         0., -1.,  2.,  0.,  1.,      -1.,    0.0,      0.,  0.0,
     /         0.,  0., -2.,  0.,  1.,      -1.,    0.0,      0.,  0.0,
     /         0.,  0.,  4., -2.,  2.,       1.,    0.0,      0.,  0.0/
      DATA XA/ 0.,  1.,  0.,  0.,  2.,       1.,    0.0,      0.,  0.0,
     /         1.,  1.,  2., -2.,  2.,       1.,    0.0,     -1.,  0.0,
     /         3.,  0.,  2., -2.,  2.,       1.,    0.0,      0.,  0.0,
     /        -2.,  0.,  2.,  2.,  2.,       1.,    0.0,     -1.,  0.0,
     /        -1.,  0.,  0.,  0.,  2.,       1.,    0.0,     -1.,  0.0,
     /         0.,  0., -2.,  2.,  1.,       1.,    0.0,      0.,  0.0,
     /         0.,  1.,  2.,  0.,  1.,       1.,    0.0,      0.,  0.0,
     /        -1.,  0.,  4.,  0.,  2.,       1.,    0.0,      0.,  0.0,
     /         2.,  1.,  0., -2.,  0.,       1.,    0.0,      0.,  0.0,
     /         2.,  0.,  0.,  2.,  0.,       1.,    0.0,      0.,  0.0/
      DATA XB/ 2.,  0.,  2., -2.,  1.,       1.,    0.0,     -1.,  0.0,
     /         2.,  0., -2.,  0.,  1.,       1.,    0.0,      0.,  0.0,
     /         1., -1.,  0., -2.,  0.,       1.,    0.0,      0.,  0.0,
     /        -1.,  0.,  0.,  1.,  1.,       1.,    0.0,      0.,  0.0,
     /        -1., -1.,  0.,  2.,  1.,       1.,    0.0,      0.,  0.0,
     /         0.,  1.,  0.,  1.,  0.,       1.,    0.0,      0.,  0.0/
C
C***********************************************************************
C
C
C     COMPUTATION OF  FUNDAMENTAL ARGUMENTS
C
C
      L = ((+0.064D0 * T + 31.310D0) * T + 715922.633D0) * T
     /     + 485866.733D0 + DMOD(1325.0D0*T,1.0D0) * 1296000.0D0
      L = DMOD(L,1296000.0D0)
C
      LP = ((-0.012D0 * T - 0.577D0) * T + 1292581.224D0) * T
     /     + 1287099.804D00 + DMOD(99.0D0*T,1.0D0) * 1296000.0D0
      LP = DMOD(LP,1296000.0D0)
C
      F = ((+0.011D0 * T - 13.257D0) * T + 295263.137D0) * T
     /     + 335778.877D0 + DMOD(1342.0D0*T,1.0D0) * 1296000.0D0
      F = DMOD(F,1296000.0D0)
C
      D = ((+0.019D0 * T - 6.891D0) * T + 1105601.328D0) * T
     /     + 1072261.307D0 + DMOD(1236.0D0*T,1.0D0) * 1296000.0D0
      D = DMOD(D,1296000.0D0)
C
      OM = ((0.008D0 * T + 7.455D0) * T - 482890.539D0) * T
     /     + 450160.280D0  - DMOD(5.0D0*T,1.0D0) * 1296000.0D0
      OM = DMOD(OM,1296000.0D0)
C
C***********************************************************************
C
C
      DPSI = 0.D0
      DEPS = 0.D0
C
C     SUM NUTATION SERIES TERMS, FROM SMALLEST TO LARGEST
C
      DO 10 J=1,106
      I = 107 - J
C     FORMATION OF MULTIPLES OF ARGUMENTS
      ARG = DBLE(X(1,I)) * L
     /    + DBLE(X(2,I)) * LP
     /    + DBLE(X(3,I)) * F
     /    + DBLE(X(4,I)) * D
     /    + DBLE(X(5,I)) * OM
      ARG = DMOD(ARG,1296000.0D0) / SECCON
C     EVALUATE NUTATION
      DPSI = (DBLE(X(6,I)) + DBLE(X(7,I))*T) * DSIN(ARG) + DPSI
      DEPS = (DBLE(X(8,I)) + DBLE(X(9,I))*T) * DCOS(ARG) + DEPS
   10 CONTINUE
C
C
      DPSI = DPSI * 1.0D-4
      DEPS = DEPS * 1.0D-4
C
C***********************************************************************
C
C
      RETURN
      END
*
*   MEMBER 'VASOL3' FOLLOWS
*
      SUBROUTINE SOLSYS (TJD,M,K,POS,VEL,IERR)
      save
C
C     SUBROUTINE SOLSYS VERSION 3.
C     THIS SUBROUTINE PROVIDES THE POSITION AND VELOCITY OF THE
C     EARTH AT EPOCH TJD BY EVALUATING A CLOSED-FORM THEORY WITHOUT
C     REFERENCE TO AN EXTERNAL FILE.  THIS ROUTINE CAN ALSO PROVIDE
C     THE POSITION AND VELOCITY OF THE SUN.
C
C          TJD  = TDB JULIAN DATE OF DESIRED EPOCH (IN)
C          M    = BODY IDENTIFICATION NUMBER (IN)
C                 SET M=0 OR M=1 OR M=10 FOR THE SUN
C                 SET M=2 OR M=3 FOR THE EARTH
C          K    = ORIGIN SELECTION CODE (IN)
C                 SET K=0 FOR ORIGIN AT SOLAR SYSTEM BARYCENTER
C                 SET K=1 FOR ORIGIN AT CENTER OF SUN
C          POS  = POSITION VECTOR, EQUATORIAL RECTANGULAR
C                 COORDINATES, REFERRED TO MEAN EQUATOR AND EQUINOX
C                 OF J2000.0, COMPONENTS IN AU (OUT)
C          VEL  = VELOCITY VECTOR, EQUATORIAL RECTANGULAR
C                 COORDINATES, REFERRED TO MEAN EQUATOR AND EQUINOX
C                 OF J2000.0, COMPONENTS IN AU/DAY (OUT)
C          IERR = ERROR INDICATOR (OUT)
C                 IERR=0 MEANS EVERYTHING OK
C
C
      DOUBLE PRECISION TJD,POS,VEL,
     /     PM,PA,PL,PN,TWOPI,TLAST,T0,OBL,SINE,COSE,TMASS,
     /     QJD,EL,C,P,
     /     F,PBARY,VBARY,DLON,SINL,COSL,X,Y,Z,XDOT,YDOT,ZDOT,
     /     DFLOAT,DABS,DMOD,DSIN,DCOS
      DIMENSION POS(3), VEL(3), EL(21), C(13), P(3,3),
     /     PM(4), PA(4), PL(4), PN(4),
     /     PBARY(3), VBARY(3)
      DATA EL,C/34*0.0D0/
C
C     ARRAYS BELOW CONTAIN DATA ON THE FOUR LARGEST PLANETS
C     THIS DATA USED FOR BARYCENTER COMPUTATIONS
C                 JUPITER        SATURN        URANUS       NEPTUNE
      DATA PM/   1047.355D 0,    3498.5D 0,    22869.D 0,    19314.D 0/
      DATA PA/   5.202803D 0,  9.538843D 0, 19.181951D 0, 30.057779D 0/
      DATA PL/   0.5999  D 0,  0.8728  D 0,  5.4714  D 0,  5.3269  D 0/
      DATA PN/   1.450216D-3,  5.839824D-4,  2.047627D-4,  1.043900D-4/
C
      DATA TWOPI/6.283185307179586D0/,     TLAST/0.0D0/
      DATA T0/2451545.000D0/,     OBL/23.43929111D0/
C     T0 = TDB JULIAN DATE OF EPOCH J2000.0
C     OBL = OBLIQUITY OF ECLIPTIC AT EPOCH J2000.0
C
C
      IF (TLAST.GT.0.0D0) GO TO 12
      SINE = DSIN (OBL * TWOPI/360.0D0)
      COSE = DCOS (OBL * TWOPI/360.0D0)
      TMASS = 1.0D0
      DO 10 I=1,4
   10 TMASS = TMASS + 1.0D0 / PM(I)
      TLAST = 1.0D0
   12 IERR = 0
      IF (TJD.LT.2340000.5D0) IERR = 1
      IF (TJD.GT.2560000.5D0) IERR = 2
      IF (IERR.NE.0) GO TO 110
      IF (M.GE.10) GO TO 20
      IF (M.GE. 2) GO TO 30
C
C     FORM HELIOCENTRIC COORDINATES OF SUN
   20 DO 25 J=1,3
      POS(J) = 0.0D0
   25 VEL(J) = 0.0D0
      IF (K) 90,90,110
C
C     FORM HELIOCENTRIC COORDINATES OF EARTH
   30 DO 35 I=1,3
      QJD = TJD + DFLOAT(I-2) * 0.10D0
C     SUBROUTINE SUN COMPUTES EARTH-SUN VECTOR
      CALL SUN (QJD,EL,C)
      CALL PRECES (QJD,C(11),T0,POS)
      P(I,1) = -POS(1)
      P(I,2) = -POS(2)
      P(I,3) = -POS(3)
   35 CONTINUE
      DO 40 J=1,3
      POS(J) =  P(2,J)
      VEL(J) = (P(3,J) - P(1,J)) / 0.20D0
   40 CONTINUE
      IF (K) 90,90,110
C
C     IF K=0, MOVE ORIGIN TO SOLAR SYSTEM BARYCENTER
C     SOLAR SYSTEM BARYCENTER COORDINATES COMPUTED FROM ROUGH
C     APPROXIMATIONS OF THE COORDINATES OF THE FOUR LARGEST PLANETS
   90 IF (DABS(TJD-TLAST).LT.1.0D-6) GO TO 99
      DO 92 J=1,3
      PBARY(J) = 0.0D0
   92 VBARY(J) = 0.0D0
C     THE FOLLOWING LOOP CYCLES ONCE FOR EACH OF THE FOUR PLANETS
      DO 98 I=1,4
      DLON = PL(I) + PN(I) * (TJD - T0)
      DLON = DMOD(DLON,TWOPI)
      SINL = DSIN(DLON)
      COSL = DCOS(DLON)
C     SINL AND COSL ARE THE SINE AND COSINE OF PLANET'S MEAN LONGITUDE
      X    =  PA(I) * COSL
      Y    =  PA(I) * SINL * COSE
      Z    =  PA(I) * SINL * SINE
      XDOT = -PA(I) * PN(I) * SINL
      YDOT =  PA(I) * PN(I) * COSL * COSE
      ZDOT =  PA(I) * PN(I) * COSL * SINE
      F = 1.0D0 / (PM(I) * TMASS)
      PBARY(1) = PBARY(1) + X * F
      PBARY(2) = PBARY(2) + Y * F
      PBARY(3) = PBARY(3) + Z * F
      VBARY(1) = VBARY(1) + XDOT * F
      VBARY(2) = VBARY(2) + YDOT * F
      VBARY(3) = VBARY(3) + ZDOT * F
   98 CONTINUE
      TLAST = TJD
   99 DO 100 J=1,3
      POS(J) = POS(J) - PBARY(J)
  100 VEL(J) = VEL(J) - VBARY(J)
C
  110 RETURN
      END
      SUBROUTINE SUN (DJ,EL,C)
      save
C
C     FOR USE WITH SUBROUTINE SOLSYS VERSION 3.
C     THIS SUBROUTINE COMPUTES THE COORDINATES OF THE EARTH-SUN
C     POSITION VECTOR WITH RESPECT TO THE ECLIPTIC AND EQUATOR
C     OF DATE.  A MODIFIED FORM OF NEWCOMB'S THEORY ('TABLES OF THE
C     SUN', 1898) IS USED.  ONLY THE LARGEST PERIODIC PERTURBATIONS
C     ARE EVALUATED, AND VAN FLANDERN'S EXPRESSIONS FOR THE FUNDAMENTAL
C     ARGUMENTS ('IMPROVED MEAN ELEMENTS FOR THE EARTH AND MOON', 1981)
C     ARE USED.  THE ABSOLUTE ACCURACY IS OF ORDER 1 ARCSECOND
C     FOR EPOCHS NEAR THE YEAR 2000.
C     (ADAPTED FROM SUBROUTINE IAUSUN BY P. M. JANICZEK, USNO.)
C
C          DJ   = TDB JULIAN DATE OF DESIRED EPOCH (IN)
C          EL   = ARRAY OF ORBITAL ELEMENTS (SEE BELOW) FOR
C                 EPOCH DJ (OUT)
C          C    = ARRAY OF COORDINATES (SEE BELOW) FOR
C                 EPOCH DJ (OUT)
C
C
      DOUBLE PRECISION DJ,EL,C,T,TP,T20,RO,GV,GM,GJ,GS,DL,DR,DB,DG,
     1 DBLARG,D,TWOPI,STR,RTD,R,TR,
     2 SINO,COSO,SINL,COSL,SINB,COSB,
     3 DSIN,DCOS,DABS,DMOD
C
      DIMENSION EL(21)
C
C     EL( 1)= SEMI-MAJOR AXIS, AU
C     EL( 2)= ORBITAL ECCENTRICITY
C     EL( 5)= LONGITUDE OF PERIGEE, RADIANS
C     EL( 9)= UNPERTURBED MEAN LONGITUDE, RADIANS
C     EL(10)= MEAN ANOMALY, AFFECTED BY LONG-PD PERTURBATIONS, RADIANS
C     EL(11)= UNPERTURBED RADIUS, AU
C     EL(12)= EQUATION OF THE CENTER, RADIANS
C     EL(13)= MEAN OBLIQUITY OF ECLIPTIC, RADIANS
C     EL(14)= MEAN LONGITUDE OF MOON, RADIANS
C     EL(15)= MEAN ANOMALY OF MOON, RADIANS
C     EL(16)= LUNAR MEAN ARGUMENT OF LATITUDE, RADIANS
C     EL(17)= MEAN LONGITUDE OF LUNAR ASCENDING NODE, RADIANS
C     EL(21)= MEAN LONGITUDE OF MOON'S PERIGEE, RADIANS
C             (REMAINING ELEMENTS OF ARRAY EL NOT USED)
C
      DIMENSION C(13)
C
C     C( 1) = PERTURBED RADIUS VECTOR, AU
C     C( 2) = SAME AS C(4), DEGREES
C     C( 3) = SAME AS C(5), DEGREES
C     C( 4) = ECLIPTIC LONGITUDE WRT MEAN ECL & EQUX OF DATE, RADIANS
C     C( 5) = ECLIPTIC LATITUDE  WRT MEAN ECL        OF DATE, RADIANS
C     C(11) = EQUATORIAL X WRT MEAN EQU & EQUX OF DATE, AU
C     C(12) = EQUATORIAL Y WRT MEAN EQU & EQUX OF DATE, AU
C     C(13) = EQUATORIAL Z WRT MEAN EQU & EQUX OF DATE, AU
C             (REMAINING ELEMENTS OF ARRAY C NOT USED)
C
C
C***********************************************************************
C
C     PART I    TABLES OF THE PERTURBATIONS
C
      DIMENSION X(8,46), X1(80), X2(80), X3(80), X4(80), X5(48)
      EQUIVALENCE (X(1, 1),X1(1))
      EQUIVALENCE (X(1,11),X2(1))
      EQUIVALENCE (X(1,21),X3(1))
      EQUIVALENCE (X(1,31),X4(1))
      EQUIVALENCE (X(1,41),X5(1))
C
C     PERTURBATIONS BY VENUS
C                  J    I     VC      VS    RHOC    RHOS      BC     BS
      DATA X1 /  - 1.,  0., +  33.,-  67., -  85.,-  39., +  24.,-  17.,
     2           - 1.,+ 1., +2353.,-4228., -2062.,-1146., -   4.,+   3.,
     3           - 1.,+ 2., -  65.,-  34., +  68.,-  14., +   6.,-  92.,
     4           - 2.,+ 1., -  99.,+  60., +  84.,+ 136., +  23.,-   3.,
     5           - 2.,+ 2., -4702.,+2903., +3593.,+5822., +  10.,-   6.,
     6           - 2.,+ 3., +1795.,-1737., - 596.,- 632., +  37.,-  56.,
     7           - 3.,+ 3., - 666.,+  27., +  44.,+1044., +   8.,+   1.,
     8           - 3.,+ 4., +1508.,- 397., - 381.,-1448., + 185.,- 100.,
     9           - 3.,+ 5., + 763.,- 684., + 126.,+ 148., +   6.,-   3.,
     *           - 4.,+ 4., - 188.,-  93., - 166.,+ 337.,     0.,    0./
      DATA X2 /  - 4.,+ 5., - 139.,-  38., -  51.,+ 189., -  31.,-   1.,
     2           - 4.,+ 6., + 146.,-  42., -  25.,-  91., +  12.,    0.,
     3           - 5.,+ 5., -  47.,-  69., - 134.,+  93.,     0.,    0.,
     4           - 5.,+ 7., - 119.,-  33., -  37.,+ 136., -  18.,-   6.,
     5           - 5.,+ 8., + 154.,    0.,     0.,-  26.,     0.,    0.,
     6           - 6.,+ 6., -   4.,-  38., -  80.,+   8.,     0.,    0.,
C
C     PERTURBATIONS BY MARS
C                  J    I     VC      VS    RHOC    RHOS      BC     BS
     7           + 1.,- 1., - 216.,- 167., -  92.,+ 119.,     0.,    0.,
     8           + 2.,- 2., +1963.,- 567., - 573.,-1976.,     0.,-   8.,
     9           + 2.,- 1., -1659.,- 617., +  64.,- 137.,     0.,    0.,
     *           + 3.,- 3., +  53.,- 118., - 154.,-  67.,     0.,    0./
      DATA X3 /  + 3.,- 2., + 396.,- 153., -  77.,- 201.,     0.,    0.,
     2           + 4.,- 3., - 131.,+ 483., + 461.,+ 125., +   7.,+   1.,
     3           + 4.,- 2., + 526.,- 256., +  43.,+  96.,     0.,    0.,
     4           + 5.,- 4., +  49.,+  69., +  87.,-  62.,     0.,    0.,
     5           + 5.,- 3., -  38.,+ 200., +  87.,+  17.,     0.,    0.,
     6           + 6.,- 4., - 104.,- 113., - 102.,+  94.,     0.,    0.,
     7           + 6.,- 3., -  11.,+ 100., -  27.,-   4.,     0.,    0.,
     8           + 7.,- 4., -  78.,-  72., -  26.,+  28.,     0.,    0.,
     9           + 9.,- 5., +  60.,-  15., -   4.,-  17.,     0.,    0.,
     *           +15.,- 8., + 200.,-  30., -   1.,-   6.,     0.,    0./
C
C     PERTURBATIONS BY JUPITER
C                  J    I     VC      VS    RHOC    RHOS      BC     BS
      DATA X4 /  + 1.,- 2., - 155.,-  52., -  78.,+ 193., +   7.,    0.,
     2           + 1.,- 1., -7208.,+  59., +  56.,+7067., -   1.,+  17.,
     3           + 1.,  0., - 307.,-2582., + 227.,-  89., +  16.,    0.,
     4           + 1.,+ 1., +   8.,-  73., +  79.,+   9., +   1.,+  23.,
     5           + 2.,- 3., +  11.,+  68., + 102.,-  17.,     0.,    0.,
     6           + 2.,- 2., + 136.,+2728., +4021.,- 203.,     0.,    0.,
     7           + 2.,- 1., - 537.,+1518., +1376.,+ 486., +  13.,+ 166.,
     8           + 3.,- 3., - 162.,+  27., +  43.,+ 278.,     0.,    0.,
     9           + 3.,- 2., +  71.,+ 551., + 796.,- 104., +   6.,-   1.,
     *           + 3.,- 1., -  31.,+ 208., + 172.,+  26., +   1.,+  18./
      DATA X5 /  + 4.,- 3., -  43.,+   9., +  13.,+  73.,     0.,    0.,
     2           + 4.,- 2., +  17.,+  78., + 110.,-  24.,     0.,    0.,
C
C     PERTURBATIONS BY SATURN
C                  J    I     VC      VS    RHOC    RHOS      BC     BS
     3           + 1.,- 1., -  77.,+ 412., + 422.,+  79., +   1.,+   6.,
     4           + 1.,  0., -   3.,- 320., +   8.,-   1.,     0.,    0.,
     5           + 2.,- 2., +  38.,- 101., - 152.,-  57.,     0.,    0.,
     6           + 2.,- 1., +  45.,- 103., - 103.,-  44.,     0.,    0./
C
C
C***********************************************************************
C
C     PART II   NECESSARY PRELIMINARIES
C
      DATA TWOPI /6.283185307179586D0/
      DATA STR   /206264806.2470964D0/
      DATA RTD   /57.295779513082321D0/
      DATA R     /1296000.0D0/
      TR = 1000.0D0 / STR
C
C     T  = TIME IN JULIAN CENTURIES FROM 1900 JANUARY 0
      T  = (DJ - 2415020.D0)/36525.D0
C
C     TP = TIME IN JULIAN YEARS     FROM 1850 JANUARY 0
      TP = (DJ - 2396758.D0)/365.25D0
C
C     T20= TIME IN JULIAN CENTURIES FROM J2000.0
      T20= (DJ - 2451545.D0)/36525.D0
C
C
C***********************************************************************
C
C     PART III  COMPUTATION OF ELLIPTIC ELEMENTS AND SECULAR TERMS
C
C     VAN FLANDERN'S EXPRESSIONS FOR MEAN ELEMENTS
      EL( 1) = 1.00000030007166D0
      EL( 2) = 0.016708320D0 + (-0.42229D-04 - 0.126D-06 * T20) * T20
      EL( 5) = 1018578.046D0 + (6190.046D0 +
     1                (1.666D0 + 0.012D0 * T20) * T20) * T20
      EL( 5) = EL( 5) * TR
      EL( 9) = 1009677.850D0 + (100.0D0 * R + 2771.27D0 +
     1                1.089D0 * T20) * T20
      EL( 9) = DMOD (EL( 9) * TR, TWOPI)
      EL(10) = 1287099.804D0 + (99.0D0 * R + 1292581.224D0 +
     1                (-0.577D0 - 0.012D0 * T20) * T20) * T20
      EL(10) = DMOD (EL(10) * TR, TWOPI)
      EL(13) = 84381.448D0 + (-46.8150D0 +
     1               (-0.00059D0 + 0.001813D0 * T20) * T20) * T20
      EL(13) = EL(13) * TR
C
C
C***********************************************************************
C
C     PART IV   LUNAR TERMS
C
C     VAN FLANDERN'S EXPRESSIONS FOR MEAN ELEMENTS
      EL(14) = 785939.157D0 + (1336.0D0 * R + 1108372.598D0
     1                + (-5.802D0 + 0.019D0 * T20) * T20) * T20
      EL(14) = DMOD (EL(14) * TR, TWOPI)
      EL(17) = 450160.280D0 + (-5.0D0 * R - 482890.539D0 +
     1                (7.455D0 + 0.008D0 * T20) * T20) * T20
      EL(17) = DMOD (EL(17) * TR, TWOPI)
      EL(21) = 300072.424D0 + (11.0D0 * R + 392449.965D0 +
     1                (-37.112D0 - 0.045D0 * T20) * T20) * T20
      EL(21) = DMOD (EL(21) * TR, TWOPI)
C
C     DERIVED ARGUMENTS
      EL(15) = EL(14) - EL(21)
      EL(16) = EL(14) - EL(17)
      EL(15) = DMOD (EL(15),TWOPI)
      EL(16) = DMOD (EL(16),TWOPI)
C     MEAN ELONGATION
      D      = EL(14) - EL(9)
C
C     COMBINATIONS OF ARGUMENTS AND THE PERTURBATIONS
      D = DMOD (D,TWOPI)
      ARG = D
      DL =    +  6469.*SIN(ARG) +  13.*SIN(3.*ARG)
      DR =    + 13390.*COS(ARG) +  30.*COS(3.*ARG)
C
      DBLARG = D + EL(15)
      DBLARG = DMOD (DBLARG,TWOPI)
      ARG = DBLARG
      DL = DL +  177.*SIN(ARG)
      DR = DR +  370.*COS(ARG)
C
      DBLARG = D - EL(15)
      DBLARG = DMOD (DBLARG,TWOPI)
      ARG = DBLARG
      DL = DL -  424.*SIN(ARG)
      DR = DR - 1330.*COS(ARG)
C
      DBLARG = 3.D0*D - EL(15)
      DBLARG = DMOD (DBLARG,TWOPI)
      ARG = DBLARG
      DL = DL +   39.*SIN(ARG)
      DR = DR +   80.*COS(ARG)
C
      DBLARG = D + EL(10)
      DBLARG = DMOD (DBLARG,TWOPI)
      ARG = DBLARG
      DL = DL -   64.*SIN(ARG)
      DR = DR -  140.*COS(ARG)
C
      DBLARG = D - EL(10)
      DBLARG = DMOD (DBLARG,TWOPI)
      ARG = DBLARG
      DL = DL +  172.*SIN(ARG)
      DR = DR +  360.*COS(ARG)
C
      EL(16) = DMOD (EL(16),TWOPI)
      ARG = EL(16)
      DB =    + 576.*SIN(ARG)
C
C
C***********************************************************************
C
C     PART V    COMPUTATION OF PERIODIC PERTURBATIONS
C
C     THE PERTURBING MEAN ANOMALIES
C
      GV  = 0.19984020D+01 + .1021322923D+02*TP
      GM  = 0.19173489D+01 + .3340556174D+01*TP
      GJ  = 0.25836283D+01 + .5296346478D+00*TP
      GS  = 0.49692316D+01 + .2132432808D+00*TP
      GV  = DMOD (GV,TWOPI)
      GM  = DMOD (GM,TWOPI)
      GJ  = DMOD (GJ,TWOPI)
      GS  = DMOD (GS,TWOPI)
C
C
C     MODIFICATION OF FUNDAMENTAL ARGUMENTS
C
C     APPLICATION OF THE JUPITER-SATURN GREAT INEQUALITY
C     TO JUPITER'S MEAN ANOMALY
C
      GJ = GJ + 0.579904067D-02 * DSIN (5.D0*GS - 2.D0*GJ
     1                 + 1.1719644977D0 - 0.397401726D-03*TP)
      GJ = DMOD (GJ,TWOPI)
C
C     LONG PERIOD PERTURBATIONS OF MEAN ANOMALY
C
      ST = T
C                ARGUMENT IS ( 4 MARS - 7 EARTH + 3 VENUS )
      DG = 266.* SIN (0.555015 + 2.076942*ST)
C                ARGUMENT IS ( 3 JUPITER - 8 MARS + 4 EARTH )
     1    + 6400.* SIN (4.035027 + 0.3525565*ST)
C                ARGUMENT IS ( 13 EARTH - 8 VENUS )
     2    + (1882.-16.*ST) * SIN (0.9990265 + 2.622706*ST)
C
C
C     COMPUTATION OF THE EQUATION OF THE CENTER
C
C     FORM PERTURBED MEAN ANOMALY
      EL(10) = DG/STR + EL(10)
      EL(10) = DMOD (EL(10),TWOPI)
      EL(12) =   DSIN(     EL(10)) * (6910057.D0 -(17240.D0+52.D0*T)*T)
     1         + DSIN(2.D0*EL(10)) * (  72338.D0 -    361.D0*T)
     2         + DSIN(3.D0*EL(10)) * (   1054.D0 -      1.D0*T)
C
C     THE UNPERTURBED RADIUS VECTOR
      RO     =                          30570.D0 -    150.D0*T
     1         - DCOS(     EL(10)) * (7274120.D0 - (18140.D0+50.D0*T)*T)
     2         - DCOS(2.D0*EL(10)) * (  91380.D0 -    460.D0*T)
     3         - DCOS(3.D0*EL(10)) * (   1450.D0 -     10.D0*T)
      EL(11) = 10.D0**(RO*1.D-09)
C
C
C     SELECTED PLANETARY PERTURBATIONS FROM NEWCOMB'S THEORY FOLLOW
C
C     PERTURBATIONS BY VENUS
      DO 20 K=1,16
C     ARGUMENT J * VENUS +   I * EARTH
      DBLARG = X(1,K)*GV + X(2,K)*EL(10)
      DBLARG = DMOD (DBLARG,TWOPI)
      ARG = DBLARG
      CS  = COS(ARG)
      SS  = SIN(ARG)
      DL  =(X(3,K)*CS  + X(4,K)*SS )+ DL
      DR  =(X(5,K)*CS  + X(6,K)*SS )+ DR
      DB  =(X(7,K)*CS  + X(8,K)*SS )+ DB
   20 CONTINUE
C
C     PERTURBATIONS BY MARS
      DO 30 K=17,30
C     ARGUMENT  J * MARS +   I * EARTH
      DBLARG = X(1,K)*GM + X(2,K)*EL(10)
      DBLARG = DMOD (DBLARG,TWOPI)
      ARG = DBLARG
      CS  = COS(ARG)
      SS  = SIN(ARG)
      DL  =(X(3,K)*CS  + X(4,K)*SS )+ DL
      DR  =(X(5,K)*CS  + X(6,K)*SS )+ DR
      DB  =(X(7,K)*CS  + X(8,K)*SS )+ DB
   30 CONTINUE
C
C     PERTURBATIONS BY JUPITER
      DO 40 K=31,42
C     ARGUMENT J*JUPITER +   I * EARTH
      DBLARG = X(1,K)*GJ + X(2,K)*EL(10)
      DBLARG = DMOD (DBLARG,TWOPI)
      ARG = DBLARG
      CS  = COS(ARG)
      SS  = SIN(ARG)
      DL  =(X(3,K)*CS  + X(4,K)*SS )+ DL
      DR  =(X(5,K)*CS  + X(6,K)*SS )+ DR
      DB  =(X(7,K)*CS  + X(8,K)*SS )+ DB
   40 CONTINUE
C
C     PERTURBATIONS BY SATURN
      DO 50 K=43,46
C     ARGUMENT J*SATURN  +   I * EARTH
      DBLARG = X(1,K)*GS + X(2,K)*EL(10)
      DBLARG = DMOD (DBLARG,TWOPI)
      ARG = DBLARG
      CS  = COS(ARG)
      SS  = SIN(ARG)
      DL  =(X(3,K)*CS  + X(4,K)*SS )+ DL
      DR  =(X(5,K)*CS  + X(6,K)*SS )+ DR
      DB  =(X(7,K)*CS  + X(8,K)*SS )+ DB
   50 CONTINUE
C
C
C***********************************************************************
C
C     PART VI   COMPUTATION OF ECLIPTIC AND EQUATORIAL COORDINATES
C
      C(1) = EL(11)*10.D0**(DR*1.D-09)
      C(4) = (DL + DG + EL(12))/STR + EL(9)
      C(4) = DMOD (C(4),TWOPI)
      C(5) = DB/STR
      C(2) = C(4)*RTD
      C(3) = C(5)*RTD
      SINO = DSIN (EL(13))
      COSO = DCOS (EL(13))
      SINL = DSIN (C(4))
      COSL = DCOS (C(4))
      SINB = DSIN (C(5))
      COSB = DCOS (C(5))
      C(11) = C(1) * (COSB * COSL)
      C(12) = C(1) * (COSB * SINL * COSO - SINB * SINO)
      C(13) = C(1) * (COSB * SINL * SINO + SINB * COSO)
C
C
C***********************************************************************
C
C
      RETURN
      END
